NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


THESIS 


A  3D  PARABOLIC  EQUATION  (PE)  BASED  TECHNIQUE 
FOR  PREDICTING  PROPAGATION  PATH  LOSS  IN  AN 

URBAN  AREA 

by 

Keem  Boon  Thiem 
September  2001 

Thesis  Supervisor:  Ramakrishna  Janaswamy 

Thesis  Committee  Members:  David  Jenn 

TriHa 


Approved  for  public  release;  distribution  is  unlimited 


Report  Documentation  Page 


Report  Date 

Report  Type 

Dates  Covered  (from...  to) 

30  Sep  2001 

N/A 

- 

Title  and  Subtitle 

Contract  Number 

A  3D  Parabolic  Equation  (PE)  Based  Technique  for 
Predicting  Propagation  Path  Loss  in  an  Urban  Area 

Grant  Number 

Program  Element  Number 

Author(s) 

Project  Number 

Keem  Boon  Thiem 

Task  Number 

Work  Unit  Number 

Performing  Organization  Name(s)  and  Address(es) 

Research  Office  Naval  Postgraduate  School  Monterey, 

Ca  93943-5138 

Performing  Organization  Report  Number 

Sponsoring/Monitoring  Agency  Name(s)  and 

Sponsor/Monitor’s  Acronym(s) 

Address(es) 

Sponsor/Monitor’s  Report  Number(s) 

Distribution/ Availability  Statement 

Approved  for  public  release,  distribution  unlimited 

Supplementary  Notes 

Abstract 

Subject  Terms 

Report  Classification 

unclassified 

Classification  of  this  page 

unclassified 

Classification  of  Abstract 

unclassified 

Limitation  of  Abstract 

UU 

Number  of  Pages 

108 

THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


REPORT  DOCUMENTATION  PAGE 


FonnApproveclOMBN(x  0704-018{i 
Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including 
the  time  for  reviewing  instruction,  searching  existing  data  sources,  gathering  and  maintaining  the  data  needed,  and 
completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any 
other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington 
headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite 
1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project 

(0704-0188)  Washington  DC  20503. _ 

2.  REPORT  DATE 

September  2001 


6.  author(S)  Keem  Boon  Thiem 


11.  SUPPLEMENTARY  NOTES 

The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official  policy  or  position  of  the 
Department  of  Defense  or  the  U.S.  Government. 


13.  ABSTRACT  (maximum  200  words)  A  mobile  radio  environment  places  fundamental 
limitations  on  the  performance  of  wireless  communication  systems.  Most  models  developed  to 
predict  propagation  path  loss  have  been  historically  performed  in  a  statistical  approach.  These 
models  are  expensive  to  develop  and  do  not  offer  the  accuracy,  computational  advantages,  and 
sufficiency  as  the  parabolic  equation  (PE).  The  goal  of  this  thesis  is  to  develop  a  3D  model 
based  on  PE  for  predicting  propagation  path  loss  in  urban  areas  on  flat  and  hilly  terrains.  The 
PE  method  offers  the  computational  advantages,  where  one  can  approximate  the  elliptic 
operator  governing  the  true  wave  behavior  by  a  much  simpler  parabolic  operator  that  permits 
marching  in  range.  Moreover  those  all-important  aspects  of  propagation  such  as  reflection  and 
diffraction  are  included  automatically  in  the  formulation.  Four  test  problems  on  flat  terrain  and 
two  test  problems  on  hilly  terrain  will  be  simulated.  For  the  flat  terrain,  the  3D  PE  model 
results  will  be  compared  with  the  two-ray,  the  four-ray,  the  UTD,  and  the  numerical  integration 
technique  results.  For  the  hilly  terrain,  the  results  of  the  3D  PE  model  will  be  compared  with 
the  UTD  and  the  numerical  integration  technique  results. 


16.  PRICE  CODE 


NSN  7540-0 1  -280-5500  Standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Std.  239-18 


20.  LIMITATION 
OF  ABSTRACT 

UL 


15.  NUMBER  OF 
PAGES 

106 


14.  SUBJECT  TERMS 

Radiowave  Propagation,  Parabolic  Equation  Techniques,  Wireless 
Communications,  Large-Scale  Propagation  Path  Loss 

18.  SECURITY 
CLASSIFICATION  OF  THIS 
PAGE 

Unclassified 


19.  SECURITY 
CLASSIFICATION  OF 
ABSTRACT 

Unclassified 


17.  SECURITY 
CLASSIFICATION  OF 
REPORT 

Unclassified 


12b.  DISTRIBUTION  CODE 


12a.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Naval  Postgraduate  School 
Monterey,  CA  93943-5000 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

N/A 


5.  FUNDING  NUMBERS 


1.  PERFORMING  ORGANIZATION 
1EPORT  NUMBER 


10.  SPONSORING  /  MONITORING 
AGENCY  REPORT  NUMBER 


4.  TITLE  AND  SUBTITLE: 

A  3D  Parabolic  Equation  (PE)  Based  Technique  for  Predicting  Propagation  Path 
Loss  in  an  Urban  Area 


3.  REPORT  TYPE  AND  DATES  COVERED 

Engineer’s  Thesis 


1.  AGENCY  USE  ONLY  (Leave  blank) 


1 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


11 


Approved  for  public  release:  distribution  is  unlimited 


A  3D  PARABOLIC  EQUATION  BASED  TECHNIQUE  FOR  PREDICTING 
PROPAGATION  PATH  LOSS  IN  AN  URBAN  AREA 

By 

Keem  B.  Thiem 

Electronics  Engineer,  GS-13,  United  States  Air  Force  Civil  Service 
B.S.E.E.,  California  State  Polytechnic  University,  Pomona,  1990 
M.S.E.E,  Naval  Postgraduate  School,  Monterey,  1993 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


ELECTRICAL  ENGINEER 


from  the 


NAVAL  POSTGRADUATE  SCHOOL 
September  2001 


Author: 


(JZaaA  )  @  • 

sy  Keem  B.  Thiem 


Approved  by: 


Ramakrisnna  Janaswamy,  Thesis  Supervisor 


c- 


David  Jenn,  Thefts  Committee  Member 


Department  of  Electrical  and  Computer  Engineering 


in 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


IV 


ABSTRACT 


A  mobile  radio  environment  places  fundamental  limitations  on  the  performance  of 
wireless  communication  systems.  Most  models  developed  to  predict  propagation  path 
loss  have  been  historically  performed  in  a  statistical  approach.  These  models  are 
expensive  to  develop  and  do  not  offer  the  accuracy,  computational  advantages,  and 
sufficiency  as  the  parabolic  equation  (PE).  The  goal  of  this  thesis  is  to  develop  a  3D 
model  based  on  PE  for  predicting  propagation  path  loss  in  urban  areas  on  flat  and  hilly 
terrains.  The  PE  method  offers  the  computational  advantages,  where  one  can 
approximate  the  elliptic  operator  governing  the  true  wave  behavior  by  a  much  simpler 
parabolic  operator  that  permits  marching  in  range.  Moreover  those  all-important  aspects 
of  propagation  such  as  reflection  and  diffraction  are  included  automatically  in  the 
formulation.  Four  test  problems  on  flat  terrain  and  two  test  problems  on  hilly  terrain  will 
be  simulated.  For  the  flat  terrain,  the  3D  PE  model  results  will  be  compared  with  the 
two-ray,  the  four-ray,  the  UTD,  and  the  numerical  integration  technique  results.  For  the 
hilly  terrain,  the  results  of  the  3D  PE  model  will  be  compared  with  the  UTD  and  the 
numerical  integration  technique  results. 
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EXECUTIVE  SUMMARY 


A  mobile  radio  environment  places  fundamental  limitations  on  the  performance  of 
wireless  communication  systems.  Most  models  developed  to  predict  the  propagation  path 
loss  have  been  historically  performed  by  statistical  approach.  These  models  are  often 
complicated  and  expensive  to  develop  and  do  not  offer  the  accuracy,  computational 
advantages,  and  sufficiency  such  as  the  parabolic  equation  (PE)  technique  or  a  ray  based 
technique. 

The  aim  of  this  thesis  is  to  develop  a  3D  model  based  on  PE  for  predicting  the 
propagation  path  loss  in  urban  areas  on  flat  and  hilly  terrains.  The  PE  method  offers  the 
computational  advantages,  where  one  can  approximate  the  elliptic  operator  governing  the 
true  wave  behavior  by  a  much  simpler  parabolic  operator  that  permits  marching  in  range. 
This  method  has  the  advantage  that  all-important  aspects  of  propagation  such  as 
reflection  and  diffraction  are  included  automatically  in  the  formulation.  Two  types  of 
terrains  are  considered,  the  flat  earth  and  the  hilly  terrain.  For  the  flat  earth  case,  four  test 
problems  are  examined.  We  compare  the  results  from  these  cases  with  the  results 
available  in  open  literature  from  the  two-ray  model ,  the  four-ray  model ,  the  UTD.  the 
theoretical  model  proposed  by  Lee,  and  the  numerical  integration  technique  presented  by 
Bertoni.  For  the  hilly  terrain  case,  two  test  problems  are  considered.  The  results  of  these 
cases  are  compared  with  the  results  presented  by  Piazzi  and  Bertoni. 

For  the  flat  earth,  the  first  test  problem  is  to  simulate  and  determine  the 
propagation  factor  F  at  the  final  range  over  flat  earth  without  obstacles  being  placed 
between  the  transmitter  and  the  receiver.  F  is  defined  as  the  normalized  field 

Hy(x,y,z)/ Hy(f.s.) |,  where  Hy(f.s.)  is  the  free- space  magnetic  field.  A 

transmitting  antenna  height  of  Ht  =  30  m,  and  the  operating  frequency  of  1  GHz  is  used. 
The  results  of  the  3D  PE  model  are  compared  with  the  two-ray  model.  This  is  the 
simplest  validation  case  we  consider. 

The  second  test  problem  is  to  simulate  and  determine  F  with  a  single  absorbing 
screen  placed  between  the  transmitter  and  the  receiver  that  has  a  height  =  50  m  and  a 

xvii 


width,  Wk  =  49.5  m.  The  screen  represents  a  building.  The  transmitting  antenna  height  is 
H,  =  60  m,  and  the  operating  frequency  is  1  GHz.  We  compare  the  results  of  the  3D  PE 
model  with  the  results  of  the  four-ray  model  which  considers  reflection  and  diffraction  in 
the  vertical  plane  joining  the  transmitter  and  receiver. 

In  the  third  test  problem,  nine  absorbing  screens  of  uniform  heights,  equal 
spacing,  variable  screen  widths  are  placed  between  the  transmitter  and  the  receiver.  A 
transmitting  antenna  height  of  H,  =  10  m,  and  operating  frequency  of  1  GHz  is  used.  This 
test  case  is  equivalent  to  the  radiowave  propagating  over  a  row  of  houses  in  an  urban 
area.  The  three  different  screen  widths  of  50  m,  25  m,  and  12.5  m  are  considered.  Each 
screen  has  a  height  of  10  m  and  is  spaced  50  m  from  the  previous  one.  Their  effects  on 
the  propagation  factor  F  are  studied.  The  propagation  factors  are  determined  at  the 
rooftops.  The  results  of  the  3D  PE  model  are  compared  with  the  UTD  results  presented 
by  Andersen  and  the  theoretical  results  proposed  by  Lee. 

The  last  test  problem  for  the  flat  earth  case  involves  120  absorbing  screens  of 
uniform  heights  and  equal  spacing  between  the  transmitter  and  the  receiver.  These 
screens  represent  row  of  buildings  or  houses.  The  transmitting  antenna  has  a  height  of  H, 
=  125  m,  and  the  operating  frequency  is  900  MHz.  The  screens  have  the  height,  =  20 
m,  and  the  width,  Wk  =  50  m.  They  are  spaced  50  m  apart.  The  propagation  factor  is 
determined  at  the  final  range.  We  compare  the  results  of  the  3D  PE  model  with  the 
results  of  the  UTD  method  and  the  numerical  integration  technique  presented  by  Bertoni. 

For  the  hilly  terrain  case,  the  first  test  problem  involves  a  single  rounded  hill  with 
multiple  absorbing  screens  of  uniform  heights,  equal  spacing  and  variable  screen  widths. 
The  three  variable  screen  widths  are  25  m,  50  m  and  100  m.  The  screen  height  Hk  =  7  m, 
and  the  spacing  of  50  m  are  maintained  constant.  Their  effects  on  the  propagation  factor 
F  are  studied.  The  transmitting  antenna  height  is  H,  =  57  m,  and  the  operating  frequency 
is  900  MHz.  The  field  strengths  are  determined  on  the  rooftops.  We  compare  the  results 
of  the  3D  PE  model  with  the  results  of  the  numerical  integration  technique  proposed  by 
Piazzi  and  Bertoni. 
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For  the  second  test  problem  in  this  class,  the  propagation  took  place  over  two  hills 
of  sinusoidal  shape  with  multiple  absorbing  screens  of  uniform  heights,  equal  spacing, 
and  variable  screen  widths.  The  three  variable  screen  widths  are  25  m,  50  m  and  100  m. 
A  screen  height  Hk  =  7  m,  and  spacing  of  50  m  is  used.  Their  effects  on  the  propagation 
factor  F  are  studied.  The  transmitting  antenna  has  a  height  of  H,  =  57  m,  and  the 
operating  frequency  is  900  MHz.  Again,  we  compare  the  results  of  the  3D  PE  model 
with  the  results  presented  by  Piazzi  and  Bertoni. 

Finally,  the  3D  model  results  based  on  PE  for  predicting  propagation  path  loss  in 
urban  areas  over  flat  and  hilly  terrains  are  simulated.  Six  different  test  problems  are 
considered.  We  compare  the  results  of  the  3D  PE  model  with  the  available  results  from 
the  literature,  and  they  show  excellent  agreement.  We  also  demonstrate  that  the  3D  PE 
model  can  support  both  flat  and  hilly  terrains  with  multiple  absorbing  screens  of  uniform 
heights,  equal  spacing,  and  variable  screen  widths. 
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I.  INTRODUCTION 


A.  BACKGROUND 

A  typical  mobile  radio  environment  in  an  urban  area  for  wireless  communication 
systems  has  no  direct  line-of-sight  path  between  the  transmitter  and  the  receiver.  The 
environment  is  so  dynamic  and  the  path  between  the  transmitter  and  the  receiver  can  vary 
drastically  from  simple  line-of-sight  to  one  that  is  severely  obstructed  by  buildings, 
mountains,  and  trees.  Hence,  the  mobile  radio  environment  places  fundamental 
limitations  on  the  performance  of  wireless  communication  systems.  Since  the 
environment  is  extremely  random,  most  analysis  done  on  it  has  been  historically 
performed  in  a  statistical  approach,  often  based  on  measurements  made  specifically  for 
an  intended  communication  system  or  spectrum  allocation  [Ref.  1], 

Currently,  there  are  large-scale  propagation  path  loss  models  being  utilized  in 
commercial  cellular  communication  systems  in  America,  Europe,  Asia,  and  around  the 
world.  Some  of  these  popular  models  include  the  Okumura  model,  the  Hata  model,  and 
the  COST-23 1-Walfish-Ikegami  model  [Ref.  2].  As  previously  mentioned,  these  models 
were  developed  by  measurements  and  statistical  analysis  made  specifically  for  frequency 
coverage  between  150  MHz  and  2  GHz,  and  are  typically  based  on  a  certain  city 
environment  like  Tokyo  or  New  York.  These  models  provide  reasonable  approximation 
for  current  cellular  communication  applications.  However,  they  are  complicated  and 
expensive  to  develop  and  do  not  offer  the  accuracy,  computational  advantages,  and 
efficiency  of  models  such  as  the  parabolic  equation  (PE),  ray  methods,  etc. 

Previously  developed  2D  parabolic  equation  based  on  the  vertical  plane  method 

does  not  account  for  the  lateral  propagation  of  waves.  Consequently,  the  mean  path  loss 

is  overestimated  and  the  standard  deviation  of  error  tends  to  be  rather  high  [Ref.  3]. 

Therefore,  a  3D  PE  formulation  based  on  the  split-step  algorithm  was  developed  and 

demonstrated  for  forward  propagation  around  perfectly  reflecting  obstacles  located  on 

ground  that  includes  the  laterally  propagating  waves  [Ref  4].  The  3D  formulation  is 

expected  to  substantially  improve  the  accuracy  relative  to  the  2D  vertical  plane  method 

for  the  test  problems.  The  goal  of  this  thesis  is  to  develop  a  3D  model  based  on  a  scalar 

1 


parabolic  equation  for  predicting  propagation  path  loss  in  urban  areas  on  flat  earth  and 
hilly  terrains. 

B.  RESEARCH  DEFINITION 

Given  an  operating  frequency  /,  the  transmitting  antenna  height  H,  and  the  3-dB 
beam  width,  and  the  terrain  profile,  we  want  to  calculate  the  field  at  a  certain  range  from 
the  transmitting  antenna.  We  use  the  parabolic  equation  (PE)  method  to  determine  the 
field  at  the  desired  location.  Although  there  are  several  computational  techniques  for 
predicting  the  field  at  a  desired  range  from  the  transmitting  antenna  including  the  ray 
tracing  approach,  and  the  integral  equation  method  etc.,  few  offer  the  computational 
advantages  of  the  parabolic  equation  (PE)  method,  where  one  can  approximate  the 
elliptic  operator  governing  the  true  wave  behavior  by  a  much  simpler  parabolic  operator 
that  permits  marching  in  range.  The  PE  method  has  the  advantage  that  all-important 
aspects  of  propagation  such  as  reflection  and  diffraction  are  included  automatically  in  the 
formulation.  However,  the  penalty  for  employing  the  PE  method  is  that  it  neglects  back 
scattering  [Ref.  5].  This  assumption  will  not  contribute  any  significant  errors  for  the 
class  of  applications  considered  in  this  thesis  since  radiowaves  predominantly  propagate 
in  the  forward  direction. 

We  consider  a  vertically  polarized  current  source  f  (z)  that  produces  only 

vertically  polarized  electric  field  component  E-  or  equivalently  a  circumferential 
magnetic  field  component.  We  look  at  the  y  component  of  the  magnetic  field  Hy  because 
the  depolarization  is  ignored.  Thus  initial  magnetic  field  H  (0,  y,  z)  at  the  location  x  =  0 
is  set  up  by  the  current  source  on  a  flat  plane,  and  it  is  directly  related  to  the  current 

source  density  J  .  It  is  easier  to  relate  the  magnetic  field  to  the  current  source  [Ref.  3,  4]. 
For  this  reason,  we  use  the  magnetic  field  in  the  3D  PE  algorithm  instead  of  electric  field; 
however,  the  electric  field  and  magnetic  field  in  the  far- zone  are  related  via  the  medium 
impedance.  The  field  is  then  split  into  even  and  odd  parts  to  ensure  that  there  are  no 
discontinuities  at  zero  crossing  which  leads  to  erroneous  results  during  numerical 
evaluation. 
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In  the  presence  of  vertical  screens  that  are  perpendicular  to  the  preferred  axis  (x- 
axis),  we  still  ignore  the  depolarization.  These  screens  represent  a  row  of  buildings  or 
houses  in  an  urban  area.  In  this  thesis,  we  assume  that  they  are  perfectly  absorbing  which 
means  that  all  incident  fields  on  them  are  zeroed  out.  Figure  1  illustrates  how  the  3D  PE 
model  algorithm  handles  the  incident  fields  on  the  screens.  First,  we  determine  the  fields 
H  (x,y,z)  in  the  absence  of  the  screens  and  then  make  local  adjustments  to  them.  After 

making  the  adjustments,  the  field  H  (y+,  y,z)  continues  to  march  in  range  where  the 

superscript  +  indicates  the  fields  immediately  after  the  screen.  This  is  then  repeated  till 
the  desired  receiver  location  is  reached. 


Original  Field 


Null  Field 


*►  Y 


Figure  1.  Fields  Incident  on  a  Perfectly  Absorbing  Screen. 

Therefore,  this  work  will  attempt  to  design  a  more  accurate  large-scale 
propagation  path  loss  model  based  on  the  scalar  3D  parabolic  equation  (PE)  over  the 
frequency  band  of  300  MHz  to  10  GHz.  The  algorithm  will  predict  the  lateral  and 
vertical  wave  propagation  in  an  urban  area  that  is  comprised  of  vertical  buildings  with 
arbitrary  cross  section  and  perfectly  absorbing  surfaces  on  flat  and  hilly  terrains.  The  3D 
PE  algorithm  also  includes  the  log-normal  shadowing  effects.  We  will  compare  the  3D 
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PE  results  to  the  two-ray  model  [Ref.  2],  four-ray  model  [Ref.  2],  the  results  given  by 
Andersen  [Ref.  6]  obtained  from  the  uniform  theory  of  diffraction  (UTD)  method,  the 
theoretical  results  given  by  Lee  [Ref.  7],  and  the  numerical  results  obtained  by  Piazzi  and 
Bertoni  [Refs.  8-10].  The  3D  PE  algorithm  is  implemented  in  MatLab. 

In  Chapter  II  we  provide  detailed  discussions  on  the  numerical  evaluation  of  the 
3D  parabolic  equation  and  the  MatLab  implementation.  The  current  source  has  vertical 
extent  along  the  "-axis  and  has  a  Gaussian  distribution.  The  source  standard  deviation, 
07 ,  is  chosen  to  produce  the  required  transmitting  antenna  3  dB  elevation  beam  width, 
®el .  A  detailed  discussion  on  the  current  source  and  its  relationship  to  the  magnetic  field 
is  also  given  in  this  chapter. 

In  Chapter  III  and  IV,  we  discuss  our  results  and  comparison  efforts  in  detail.  We 
compare  the  3D  parabolic  equation  (PE)  results  with  the  results  reported  in  the  literature. 
In  Chapter  III,  first,  we  examine  a  flat  earth  and  perfectly  reflecting  ground.  The  results 
of  the  3D  PE  model  of  this  case  are  compared  with  the  results  of  the  two-ray  model. 
Then  we  introduce  single  and  multiple  perfectly  absorbing  screens  of  uniform  heights  and 
equal  spacing  between  the  transmitter  and  the  receiver  over  a  flat  earth  and  non-perfectly 
reflecting  ground.  These  absorbing  screens  represent  row  of  buildings  of  arbitrary 
widths.  The  effects  of  finite  screen  widths  on  simulation  are  also  examined. 

In  Chapter  IV,  we  introduce  a  rounded  hill  and  two  hills  of  sinusoidal  shape  with 
multiple  absorbing  screens  the  uniform  heights  and  equal  spacing  between  the  transmitter 
and  the  receiver.  These  scenarios  are  representative  of  well  built-up  urban  areas.  In  both 
cases,  we  also  examine  the  effects  of  finite  screen  widths  versus  the  infinitely  long  screen 
widths  on  the  simulation  results.  Results  from  these  cases  are  compared  with  the  results 
presented  by  Piazzi  and  Bertoni.  Chapter  V  provides  the  conclusions  and 
recommendations.  All  MatLab  codes  are  listed  in  the  Appendices. 
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II.  NUMERICAL  EVALUATION  OF  3D  PARABOLIC 

EQUATION 


A  detailed  method  for  evaluating  the  3D  parabolic  equation  (PE)  numerically  is 
presented  in  this  chapter.  First  the  3D  PE  model  algorithm  flow  diagram  is  presented  in 
Section  A.  Then  the  initial  current  source  and  its  characteristics  are  discussed  in  Section 
B.  This  is  followed  in  Section  C  by  a  detailed  discussion  on  how  the  magnetic  field  is 
related  to  the  current  source.  A  Hanning  window  is  used  throughout  the  algorithm  to 
contain  the  fields  both  in  spatial  and  wavenumber  domains  before  taking  2D  (Ny  x  Nz) 
Fourier  transforms  (2FFT’s)  or  2D  inverse  Fourier  transforms  (IFFT’s).  The  Hanning 
window  is  discussed  in  detail  in  Section  D.  Brief  descriptions  of  the  auxiliary  parameters 
are  provided  in  Section  E. 

In  the  parabolic  equation  the  fields  are  coupled  causally  from  one  range  to  the 
next.  There  is  no  path  for  waves  to  return  information  from  obstacles  located  ahead  of 
the  receiver.  The  initial  fields  are  set  up  by  the  current  source  on  an  impedance  plane. 
Presence  of  a  vertical  screen  perpendicular  to  the  preferred  axis  (a- ax  is)  is  handled  by 
first  determining  the  fields  in  the  absence  of  the  screen  and  then  making  local 
adjustments  to  them.  This  process  is  then  repeated  till  the  desired  receiver  location  is 
reached. 

Figure  2  illustrates  the  geometry  of  the  scalar  3D  PE  model  approach.  To  find  the 
field  over  an  impedance  plane  at  a  particular’  range,  a  +  Ax,  given  the  field  at  a  previous 
range,  x,  the  latter  is  first  decomposed  into  a  spectrum  of  plane  waves  traveling  in  the 
positive  v-direction.  Each  plane  wave  is  then  propagated  to  the  new  range  using  the 
appropriate  propagator.  The  plane  waves  at  the  new  range  are  then  added  to  compose  the 
field  at  x  +  Ax.  The  various  spectral  decompositions  are  carried  out  using  2D  Fourier 
transforms  ( IFFTs )  in  the  y-z  plane.  Furthermore,  the  magnetic  field  can  be  defined  in 
terms  of  odd  and  even  parts  to  ensure  that  there  are  no  discontinuities  at  zero  crossing 
which  leads  to  erroneous  results  during  the  2D  Fourier  transform  ( FFT)  computation.  In 
this  thesis,  the  current  source  is  assumed  to  be  along  the  positive  "-direction,  and  the 
antenna  is  vertically  polarized.  Over  flat  plane,  the  only  non-zero  field  components  are 
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Ez,E, a  and  HCh  Where  p,  z,  (pare  the  standard  cylindrical  coordinates,  and  (f) is  measured 
from  the  x-axis.  We  ignore  depolarization  and  continue  this  assumption  to  hold  in  the 
presence  of  the  screens.  Thus  the  vertically  polarized  current  source  f  (z)  produces 

only  vertically  polarized  electric  field  component  Ez  or  equivalently  a  circumferential 
magnetic  field  component.  In  this  thesis,  we  choose  to  use  the  y  component  of  the 
magnetic  field  Hy  because  the  antenna  depolarization  is  ignored. 


Figure  2.  Geometry  of  the  Scalar  3D  PE  Model  [From:  Ref.  11]. 

A.  3D  PARABOLIC  EQUATION  ALGORITHM 

This  section  provides  an  overall  view  of  the  3D  PE  model  algorithm  flow  diagram 
as  illustrated  by  Figure  3.  The  method  of  implementing  the  algorithm  in  MatLab  is  also 
presented.  Obstacles’  placements  and  terrains  are  also  discussed  in  the  context  of 
simulation.  Additional  details  are  provided  in  the  subsequent  sections. 
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Apply  Hanning  window 
in  frequency  domain 


e 


ikx Ax  (free  space 
propagator) 


2D  IFFT  w.r.t.  ky  and  k, 

Hanning  window  in  spatial  domain 


Figure  3.  3D  PE  Numerical  Evaluation  Flow  Diagram. 


The  initial  magnetic  field  H(0,y,z )  is  assumed  to  be  generated  by  the  current 


source  on  a  flat  plane.  The  Gaussian  current  source  along  the  z-axis  gives  rise  to  a 
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circumferential  magnetic  field  with  a  lateral  component,  H y[x,y,z)  ■  It  is  described  in 

Section  C.  The  transformed  magnetic  field  H  y  (0',/vv,/v.  )  provides  the  initial  field  for 

the  beginning  of  the  main  loop  of  the  algorithm  as  shown  in  Figure  3,  where  the  tilde 

indicates  the  transformed.  When  H y  (()  ,ky,kz  )  is  multiplied  by  the  free-space 

propagator  e'k'A' ,  the  initial  transformed  field  has  marched  from  the  previous  location  x 
to  a  new  location  x  +  Ax  and  is  given  by  equation  (2.1).  The  field  is  split  into  even  and 
odd  parts  for  numerical  convenience. 

Hy  (v  +  Ax,ky,kz )  =  <  >elk^Ax  (2.1) 

The  quantity  Fy  (k  )  is  the  reflection  coefficient  of  a  plane  wave  for  parallel  polarization. 

Then  the  magnetic  field  H  (x,y,z)  at  a  new  location  v  +  Ax  just  right  before  the 
screen  is  found  by  double  integrating  equation  (2.1)  and  is  expressed  by  equation  (2.2). 

Hy  (x,  y,z)  =  — *-y  J  |  Hy  (.x  +  Ax,  ky,kz ) e‘(kyZ+k,  'dkydkz  (2.2) 

If  we  examine  equation  (2.2),  it  is  basically  a  2D  inverse  Fourier  transform,  (IFFT),  of 
the  left  hand  side  of  equation  (2.1)  with  respect  to  ky  and  kz.  Then  instead  of  evaluating 
equation  (2.2)  by  complex  double  integrations,  the  2D  (Ny  x  N-)  IFFT  is  used  to  evaluate 
it.  The  equivalent  MatLab  expression  to  equation  (2.2)  is  defined  by  equation  (2.3). 

A  k  A  k  (  ~  \ 

tf,(w)  =  7rV/FfT2  HAx’krK)  PA).  (2.3) 

(2  71)  V  J 

Also,  MatLab  divides  the  results  by  Ny  x  N-  when  it  performs  a  2D  IFFT  operation; 
therefore,  the  results  are  de-normalized  by  multiplying  them  by  Ny  x  N-.  The  results  are 
multiplied  by  a  Hanning  window  in  spatial  domain  to  prevent  aliasing  effects. 
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At  this  point,  after  the  magnetic  field  H  y  ( x,y,z )  has  been  found  just  before  the 

screen,  the  obstacles  are  introduced  that  represent  buildings,  houses,  trees,  etc.  in  the 
simulation.  In  all  the  test  problems,  screens  are  used  as  obstacles  and  assumed  to  be 
perfectly  absorbing.  These  screens  represent  rows  of  buildings  having  height  Hk  and 
width  Wi...  The  Hk  and  Wk  parameters  are  used  to  study  the  effects  on  the  propagation 
factor  and  field  strengths.  All  incident  fields  on  the  screens  are  eliminated  since  they  are 
perfectly  absorbing.  The  fields  that  propagate  forward  consist  of  the  diffracted  and  the 
reflected  waves.  The  H y[x+  ,y,z)  indicates  the  field  after  adjusting  for  screens  at 

location  x.  Then  the  fields  begin  their  march  immediately  after  the  screen.  Figure  4 
illustrates  an  example  of  a  hilly  terrain  with  multiple  screens  of  uniform  heights  and 
equal  spacing  between  the  transmitter  and  the  receiver.  The  propagation  factor  F  and  the 
field  strengths  may  be  determined  at  each  marching  step  as  the  fields  propagate  forward. 
The  data  are  stored  as  an  array. 


Figure  4.  Hilly  Terrain  and  Multiple  Equal  Height  and  Equally  Spaced  Screens. 

Once  the  magnetic  field  H  ( x,y,z ),  at  a  new  location  x,  has  been  found,  it  may 

be  expressed  in  terms  of  odd  and  even  parts  to  ensure  there  are  no  discontinuities  at  zero 
crossing  that  cause  erroneous  results  during  the  2D  Fourier  transform  ( FFT)  operations. 

The  odd  H  Q  (v,  y,  z)  and  even  H  ( x,y,z )  fields  are  given  as: 
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Hyo{X^^) 


Hy(x,y,z )  z>0 
-Hy(x,y,-z)  z<  0 


Hye{X^^)  = 


Hy(x,y,z)  z>  0 
Hy(x,y,-z )  z  <  0 


The  transformed  quantities  H  yo  [x,ky,kz )  and  H  ye(x,ky,kz^  are  found  by  the 


following  equations: 


H yo  ( x,  ky  ,kz)=  £  Hy0  ( x,  y',z')e  k,y  +Kz  ^ dy'dz ' 
H ye  (x,  ky  ,kz)  =  £  Hye  ( X,y',z')e  ^  +kz<  dy ' dz ' 


Equations  (2.6)  and  (2.7)  are  essentially  the  2D  Fourier  transforms,  2D  FFT’s,  of 
H v0  ( .r,  y,z)  and H ve  (.r,  y,  z)  with  respect  to  y  and  z.  The  equivalent  Matlab  operations 
of  equations  (2.6)  and  (2.7)  are  defined  as: 

H  yo{x,ky,kz)  =  AyAzFFTl(Hy0(x,y,z))  (2.8) 

H  ye  ( -V,  ky  ,kz)  =  AyAz.FFT 2  ( Hye  (x,  y ,  z ) )  (2.9) 

Next  the  odd  and  even  transformed  magnetic  fields  are  multiplied  by  the 
Hanning  window  in  frequency  domain  to  prevent  aliasing.  Then  the  odd  transformed 

H yo[x,ky,k7^is  multiplied  by  a  factor— [j  - FM  (k.  )J  ,  and  the  even  transformed 

H  ye(x,ky,kz^  is  multiplied  by  —  [j.  +  FH  [k.  )J  .  These  fields  are  then  combined  to 
provide  the  newly  transformed  magnetic  field  immediately  after  the  screen,  i.e. 

Hy  (()  ,ky,k  | .  The  steps  in  the  main  loop  as  illustrated  in  Figure  3  are  repeated  until 
the  desired  receiver  location  is  reached.  The  propagation  factor  F  at  the  final  range  is 
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defined  as  the  normalized  field  H  (x,y,z  )/ H y(f.s.)  ,  where  Hy(f.s.)  is  the  free- 
space  magnetic  field  and  is  defined  as: 

1  -  eik°R 

Hx{f.s)  =  -ik^in(O)cos{0)—g(k()cos{O))——  (2.10) 

4/r  K 

where  #and  (t)  arc  the  angles  in  spherical  coordinates,  i  =  V— 1 ,  and  g  (k  )  defined  by 

equation  (2.20)  is  the  vertical  plane  response  of  the  current  source.  R  is  a  range  between 
the  transmitter  and  receiver,  and  ko  is  the  free  space  wave  number.  R,  p,  0,  and  (j)  are 
defined  by  the  following  equations: 


=  tJ(z-H,)2  +  p2 

(2.11) 

p  =  V*2  +  r 

(2.12) 

6  —  sin-1  f — 1 

(2.13) 

UJ 

0  =  sin  — 

(2.14) 

VP) 

The  current  density  and  its  characteristics  are  discussed  in  the  next  section. 

B.  GAUSSIAN  CURRENT  SOURCE 

The  fields  are  excited  by  a  vertically  polarized  current  source  with  a  prescribed 
aperture  distribution  f  (z) .  We  choose  a  Gaussian  distribution  current  source  to  allow  us 

to  match  the  source  standard  deviation,  <J_ ,  to  the  transmitting  antenna  3  dB  elevation 
beam  width,  ®c/  [Ref.  3,  4].  It  is  also  convenient  to  relate  the  magnetic  field  to  the 
current  source.  We  consider  a  transmitting  antenna  at  height  z  =  H,  and  the  current 

source  density  J  as  defined  by  equation  (2.15). 

— >  A 

J  =  zljS(y-O)  f  (z)  (2.15) 

where 
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1 


(2.16) 


/(z)  =  - 


-{z~H, filar 


,yf2 ~n 


is  a  Gaussian  current  source  with  a  standard  deviation  of  <7  ,  Iol  is  a  current  moment 

[Ref.  11],  and  S(y  —  0)  is  a  dirac  function.  Figure  5  shows  the  vertically  polarized 

— > 

Gaussian  current  density  function  J  placed  over  an  impedance  plane.  The  current  is 
assumed  to  be  independent  of  y  and  has  a  peak  at  z  =  Ht.  (7  is  its  standard  deviation  and 
may  be  chosen  from  Table  1  to  fit  the  required  transmitting  antenna  3  dB  elevation  beam 


width,  0(/  [Ref.  3].  It  may  be  observed  that  R  —  —  Ht)  +  p2  is  the  direct  distance 

V2  2 

x“  +  y~  is  the  horizontal  distance  from  the  source,  and  0  and  (j) 


are  angles  in  spherical  coordinate  defined  by  G  =  sin  1  —  and  (f)  =  sin  1 


r  P} 

kRj 


V| 

p J 


.  For  the 


geometry  shown,  the  only  non-zero  field  components  are  Ez ,  Ef>  and  //, 


/i 


Figure  5.  Gaussian  Density  Current  Source  over  an  Impedance  Plane. 
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Figure  6  illustrates  the  Gaussian  current  source  function  f  (z)  with  H,  =  10  m  and 


<7_  =  X.  We  discuss  the  magnetic  field  in  detail  in  the  next  section. 


®el  >  de§ 

o_  IX 

60 

0.20 

45 

0.30 

30 

0.45 

15 

1.00 

10 

1.52 

5 

3.03 

Table  1 .  Transmitting  Antenna  Elevation  Beam  Width  Versus  <7  TFrom:  Ref.  31. 


The  function  f  (z)  can  be  written  in  terms  of  its  odd  and  even  functions,  f0(z) 
and  fe  (z) ,  respectively,  to  permit  a  more  accurate  computation  of  the  Fourier  transform 
in  the  vertical  direction.  The  odd  and  even  extensions  of  f  (  z)  are  given  in  equations 


(2.17)  and  (2.18). 


\f(z)  Z>0 

[-/(-z)  z  <  0 


fe(z) 


j/(z)  Z>0 

l/(-z)  Z<0 


(2.17) 

(2.18) 


The  functions  fo  (z)  and  fe(z)  are  shown  in  Figures  7  and  8.  If  we  add  th e/a(z) 

and fe(z)  together  and  divide  by  Vi,  we  will  get  back  the  original  function  /(z). 

Hence,  the  function  decomposition  into  its  odd  and  even  parts  to  facilitate  computation  is 
a  valid  one. 


f0(z) 

Figure  7 .  Odd  Pail  of  /  ( Z ) . 


14 


We  need  to  examine  the  Fourier  transform  of  the  current  source  density  J  .  We 

shall  use  the  transformed  information  to  derive  the  initial  transformed  magnetic  field  in 

— > 

Section  B.  Since  J  is  independent  of  y,  if  we  take  the  Fourier  transform  of  equation 
(2.15)  with  respect  to  z,  which  is  tantamount  to  taking  the  Fourier  transform  (FFT)  of 
f  (z) ■  The  FFT  of  8[y  —  0)  with  respect  to  y  is  equal  to  1.  Therefore,  the  Fourier 

transform  of  f  (z)  with  respect  to  z  is  given  below: 


f(K) 


—ik.H,  - 

e  7  ' e 


■Kg!  12 


§(K)e 


-ik.H, 


(2.19) 


where  g(  k.  )  is  defined  as 


*(*,) 


r/2 


(2.20) 
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It  is  easy  to  see  that 


fe{K)  =  2g(K)COS(KH ,) 

(2.21) 

fo{K)=-2is{K)sin(KH') 

(2.22) 

C.  MAGNETIC  FIELDS 


In  this  section,  we  derive  the  initial  transformed  magnetic  field  H y  (x, kY, k, ) 


from  the  transformed  current  source  f  (k  ) .  The  fields  march  in  the  transformed  space, 

and  the  propagation  factor  or  the  field  strength  is  determined  in  the  physical  space.  Then 
the  basic  3D  parabolic  equation  (PE)  algorithm  will  start  with  the  initial  transformed 

magnetic  field  H y  (0+,ky,kz )  as  illustrated  in  Figure  3,  where  0+  indicates  the  field 
immediately  after  the  screen.  Since  there  is  no  screen  for  the  initial  field,  0+  is  equal  to  0. 
First  we  need  to  define  the  relationship  between  the  initial  magnetic  field 

Hy  (0  ,  y,  z)  and  the  current  density  function  J  .  The  initial  H  (()  ,  y,  z)  relates  to  the 


current  density  J  is  described  by  equation  (2.23)  [Ref.  11] 

tf,(0+.y.z)  =  -^  j ~\S(y'-0)f{z')dy'dz' 

-°o  0 

P'1"’1'"  (*=  ~j<tyrjf(z')dz ' 

-(F  Res  [r„  (K )] )  dy  jdz  'S(y 0)/(z  V‘>‘’ 

0  (2.23) 

=  ^f(z)  +  ^]f(z')dz'. 

Z  Z  0 

^  (k,  )dh  -;FRes[r"  (t,  )]/--'=*=■' 
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where  the  second  part  of  equation  (2.23)  describes  the  surface  waves.  In  many  long- 
range  propagation  problems,  particularly  in  the  VHF  band  and  above,  it  is  possible  to 
completely  ignore  the  surface  wave  terms  as  they  decay  rapidly  with  range  [Ref.  5]. 

Then  equation  (2.23)  is  reduced  to  a  simple  relationship  between  the  Hy  (0+,y,zj  and 

the  Gaussian  current  source  function  f(z). 

Again,  to  ensure  there  is  no  discontinuity  at  zero  crossing  that  gives  erroneous 
results  in  the  3D  PE  numerical  calculation,  the  magnetic  field  may  be  defined  in  terms  of 
its  odd  and  even  parts,  H  ( x,y,z )  and  He  (v,  y,  z) ,  respectively.  The  methods  are  the 
same  as  in  equations  (2.4)  and  (2.5). 

From  equations  (2.4),  (2.5),  (2.23),  and  the  transformed  odd  fa(kz)  and  even 
fe(kz )  current  sources  as  defined  by  equations  (2.21)  and  (2.22),  we  can  easily  derive 
the  initial  transformed  odd  H yo  (0+, ky,k_ )  and  even  H  ye  [0+,ky,k7  J  magnetic  fields. 
The  initial  odd  transformed  magnetic  field  H yo  (0 +  ,ky,k_  j  is  defined  by  equation  (2.24) 

s,o(o+,*„*;) =!/„(*,) 

V  '  ’  2  (2.24) 

=  -ie~k^n  sin  (k_Ht ) 


The  initial  even  transformed  magnetic  field  H ye  (0+,ky,k_ )  is  defined  by  equation 
(2.25) 

v  y  z,  2  (2.25) 

=  12  cos  (kzHt) 


These  fields  are  defined  as  column  vectors  having  Nz  elements.  From  Section  A,  we 

— > 

know  that  the  current  density  function  J  is  independent  of  y,  and  the 
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FFT(j(y-0))  =  l;  then  we  can  repeat  the  column  of  the  initial 

transformed  H  vo  ( 0+ ,  ky ,  k. )  and H ye  ( 0+ ,  ky ,  k  j  Ny  times.  Now  we  have  the  2D  initial 

transformed  H  vo  ( 0 ' , ky , k. )  and  H ye  \0+  ,ky,k7  )  with  dimension  of  Ny  x  N-.  The 
number  of  elements  along  the  y-axis  and  z-axis  are  Ny  and  Nz,  respectively. 


The  transformed  magnetic  field  H  y  (x,ky,k7  j  is  defined  by  equation  (2.26). 

Hy(x,ky,kz)=  ]dy']dz'e-^y'[e-i^+ril(kz)ei^']Hy(x,y\z')  (2.26) 

0 


Recall  that  e  >klZ  and  e'kzZ  are  the  free-space  propagators  for  direct  and  reflected  rays, 
respectively,  and  y’  and  z’  are  dummy  variables  of  integration.  FM  ( k, )  is  the  reflection 
coefficient  for  the  parallel  polarization  as  a  function  of  kz  in  the  complex  k-  plane  and  in 
terms  of  normalized  impedance  Zs.  Fri  (k  )  and  Zs  are  expressed  as  : 


r„(*J 


K~Kzs 
K+K  Zs 


z  = 


(2.27) 

(2.28) 


£  -  £  + 


ilScr(mS  /  m ) 
f(MHz) 


(2.29) 


where/is  the  operating  frequency  measured  in  MHz.  We  use  equations  (2.4)  and  (2.5)  to 
decompose  the  magnetic  field  Hy(x,y',z')  into  Hyo(x,y',z')  and  Hye(x,y',Z)  and 

substitute  them  into  equation  (2.26).  Now  we  add  th e//w  (x,y',z')and  H  (x,y',z') 

together  for  z  >  0  and  perform  variable  substitutions  to  change  the  inner  limit  of 
integration  from  0  to  °°  to  -°o  to  °°  and  divide  it  by  Vi,  and  similarly  for  z  <  0.  The  result 
is: 
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h  ,  (x,  krkz ) = i[i + r,  (kz  )]££»*  (*,  y  u  > 


-i(kyy'+kzz'} 


dy'dz' 


}  ~  rn  (K )]  £  £  H yo  (x,y',z  ')e 


-i(kyy'+kzz') 


dy'dz' 


(2.30) 


If  we  examine  equation  (2.30),  we  see  that  the  double  integrations  of 
Hye(x,y',z')  and  Hyo(x,y',z')  are  the  2D  Fourier  transforms  ( FFT’s )  of  these  fields. 
Thus  equation  (2.30)  can  be  evaluated  with  the  2D  (Ny  x  Nz)  FFT’s  instead  of  complex 

double  integrations.  Thus  FI y  (x,ky,k7  )can  be  defined  in  terms  of  H yo(x,y,z)  and 


Hye[x,  y,z)  •  This  is  how  we  are  going  use  the  2D  FFT’s  in  the  MatLab  codes  to 
implement  equation  (2.30). 

After  making  the  adjustment  for  the  screens,  the  fields  begin  their  march  in  range 
immediately  after  the  screen  expressed  as 


Hy(x+,ky,kz 


)=|[i+r ,(K)]h  „{*\ky,kt) 


(2.31) 


Equation  (2.31)  provides  the  initial  transformed  magnetic  field  for  the  beginning  of  the 
main  loop  of  the  3D  PE  algorithm  as  shown  in  Figure  3.  Once  again,  v+  indicates  the 

fields  immediately  after  the  screen.  The  initial  fields  H yo  (x+,ky, k, )  and 

H  ye  ( A'+ ,  ky ,  k. )  are  defined  as  matrices  having  Ny  x  N-  elements .  These  fields  will  allow 

the  algorithm  to  account  for  all  the  vertical  and  lateral  propagation  of  waves,  which  is  the 
fundamental  goal  of  the  3D  PE  formulation.  Previously  developed  2D  PE  algorithms 
based  on  the  vertical  plane  method  did  not  account  for  the  laterally  propagating  of  waves. 
Consequently,  the  mean  path  loss  was  overestimated  and  the  standard  deviation  of  error 
tends  to  be  higher  [Ref.  11].  Therefore,  the  3D  PE  formulation  is  expected  to 
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substantially  improve  the  accuracy  relative  to  the  2D  vertical  plane  method  for  the  test 
problems. 

When  we  previously  performed  column-repeating  the  initial  H  vo  (0+,ky,k_ ) 
andH ye(0+  ,ky,kz )  by  Ny  times  to  create  the  2D  (Ny  x  Nz)  H  yo  (()  ,ky,k7  )  and  2D 

H ye  (0+,ky,kz ) ,  it  was  to  ensure  that  the  matrices  always  have  dimension  of  Ny  x  N- 

since  this  is  a  general  aperture  case.  When  Ny  is  not  equal  to  Nz,  we  have  a  rectangular 
aperture,  and  when  they  are  equal  to  each  other,  we  have  a  square  aperture.  The  3D  PE 
algorithm  is  implemented  to  support  any  Ny  x  Nz  aperture.  Where  Ny  and  N-  are  elements 
along  the  vertical  and  horizontal  ranges,  respectively,  with  a  suitable  power  of  2,  i.e.,  Ny 
=  1024  or  2048. 

Discussion  thus  far  pertains  to  what  happens  outside  of  the  main  loop  of  the 
algorithm.  Equation  (2.31)  provides  the  initial  transformed  field  for  the  beginning  of  the 

main  loop  for  the  3D  PE  algorithm;  then  we  multiply  H  (0+,ky,k7  )  by  the  Hanning 

window  in  frequency  domain  before  marching.  We  always  multiply  the  fields  by  the 
Hanning  window  before  taking  2D  Fourier  transform  (FFT’s)  or  2D  inverse  Fourier 
transform  (IFFT’s)  to  contain  them  in  space  and  wavenumber  domains.  The  Hanning 
window  is  discussed  in  the  next  section. 

D.  THE  HANNING  WINDOW 

In  order  to  contain  the  fields  in  space  and  wavenumber  domains  we  multiply  them 
with  the  Hanning  window  [Ref.  3]  before  taking  2D  Fourier  transform  (FFT’s)  or  2D 
inverse  Fourier  transform  (IFFT’s).  The  Hanning  window  provides  a  gradual  rolloff  to 
zero  over  the  last  quarter  of  the  domain,  and  may  be  constructed  from  the  two  Hanning 
sequences  below.  Equations  (2.32)  and  (2.33)  are  the  Hanning  sequences  along  the 
horizontal  range  ^’-direction)  and  the  vertical  range  (z-direction),  respectively.  The 
mirror  images  of  hy  (?)  and  h  (?)  below  ?  =  0  are  used  for  negative  wave  numbers. 
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for  0  <t<3Ny/8 


n 


MO 


sin2 


f  \  M 


v  Ny  J 


for  3NJS<t<Nv/2 


K  (0 


sin’ 


KNz  J 


for  0 <t<3N./S 
for  3NJS<t<Nz/2. 


(2.32) 


(2.33) 


The  Hanning  window,  in  general,  has  a  rectangular  shape  since  Ny  and  Nz  are  not 
necessarily  equal.  When  Ny  and  N-  are  equal,  then  the  Hanning  window  has  a  square 
shape  as  discussed  in  the  previous  section.  Figures  9  and  10  illustrate  the  examples  of 
512  elements  of  hx  (?)  and  h,  (?)  Hanning  sequences. 


Figure  9.  hy  (?)  Hanning  Sequence. 
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500 


Figure  10.  h,  (?)  Hanning  Sequence. 

As  in  the  previous  section,  the  Hanning  window  must  have  a  dimension  of  Ny  x 
Nz.  Again,  the  procedure  is  similar  to  Section  C;  we  first  create  a  matrix  from  the 
Hanning  sequence  h_  (?)  by  column  repeating  the  sequence  Ny  times.  This  method 

creates  a  Ny  x  Nz  matrix  from  the  h.  (?)  sequence.  Then  we  construct  the  second  matrix 
from  the  Hanning  sequence  the  h  (?)  by  row-repeating  it  Nz  times  (thereby  constructing 
a  Ny  x  Nz  matrix  from  hy(t )  sequence).  We  column-repeated  h  (?)  Ny  times  to  also 
ensure  that  we  do  not  create  a  matrix  that  has  a  column  dimension  greater  than  Ny.  The 
same  is  true  for  the  h  (?)  case.  We  row-repeated  N-  times  to  ensure  that  the  row  of  the 

matrix  does  not  have  dimension  greater  than  Nz.  This  implements  the  Hanning  window 
in  MatLab. 

Figure  11  shows  the  two  Hanning  matrices  created  from  /?-(?)  and  hv  (?) .  To 

construct  the  final  Hanning  window,  the  two  matrices  are  multiplied  together  element  by 

element.  Figure  12  shows  an  example  of  a  512  x  512  Hanning  window.  A  similar 
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Hanning  window  is  used  in  the  3D  PE  algorithm  depending  on  the  dimensions  of  Ny  and 
Nz. 


Figure  11.  h  (?)  and  li.  (t)  Hanning  Matrices. 
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Figure  12.  3D  Hanning  Window. 
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E.  SPATIAL  AND  FREQUENCY  PARAMETERS 

The  values  of  Ax,  Ay,  Az,  Akx  and  Akz,  and  the  range  of  kx,  ky,  and  k-  are  dictated 
by  the  Nyquist  sampling  theorem.  Ax  is  the  range  increment,  Ay  is  the  horizontal 
increment,  and  Az  is  the  vertical  increment.  kx ,  ky,  and  k-  are  the  range  wavenumber,  the 
horizontal  wavenumber,  and  the  vertical  height  wavenumber,  respectively.  Aky  and  Ak- 
are  the  increments  of  the  horizontal  and  the  vertical  wavenumbers,  respectively. 
References  3,  11,  and  12  provide  details  to  obtain  the  values  for  these  parameters. 

For  all  the  test  problems  in  this  thesis,  we  assume  that  Ay  =  Az  =  A,  and  Gz  =  X. 
We  pick  Ax  between  25  m  and  150  m.  If  we  choose  a  smaller  Ax,  the  computation  takes 
longer.  If  we  choose  a  larger  Ax,  the  computation  is  accelerated,  but  we  might  miss  the 
obstacles  that  are  within  in  the  marching  step.  Therefore,  Ax  has  to  be  optimally  chosen 
accordingly  for  each  problem.  However,  Ax  can  be  chosen  to  have  variable  values  within 
a  simulation.  Furthermore,  for  oz  =  X  corresponding  to  the  transmitting  antenna  of  a  3- 
dB  bandwidth  of  15°  is  shown  in  Table  2.  The  value  of  oz  can  be  picked  from  Table  2  to 
match  the  transmitting  antenna  3-dB  bandwidth. 

We  compute  the  propagation  factor  F  at  the  final  range  as  well  as  at  each 
marching  step.  We  also  consider  the  field  strengths  at  each  marching  step  and  on  the 
rooftops.  In  Chapter  El,  we  shall  present  the  results  for  the  flat  earth  case.  The  hilly 
terrain  test  problems  are  presented  in  Chapter  IV. 


25 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


26 


III.  FLAT  EARTH  RESULTS 


In  this  chapter,  we  present  results  for  the  propagation  factors  for  the  flat  earth  and 
perfectly  reflecting  ground.  The  propagation  factor  at  the  final  range  without  obstacles 
between  the  transmitter  and  the  receiver  is  simulated  in  Section  A.  The  result  of  the  3D 
PE  model  is  compared  with  the  result  of  the  two-ray  model  [Ref  2].  Next,  in  Section  B, 
we  place  a  single  absorbing  screen  between  the  transmitter  and  the  receiver  and  present 
the  result  for  the  propagation  factor  at  the  final  range.  The  screen  represents  a  building. 
The  simulation  result  is  compared  with  the  result  of  the  four-ray  model  [Ref.  2].  Then, 
in  Section  C,  we  place  nine  absorbing  screens  of  uniform  heights,  equal  spacing,  and 
variable  widths  between  the  transmitter  and  the  receiver.  These  screens  represent  a  row 
of  buildings  or  houses  in  residential  areas  of  a  city.  In  this  case,  the  relative  propagation 
factors  are  measured  at  the  rooftops  for  each  marching  step.  The  results  of  the  3D  PE 
model  are  compared  with  the  results  presented  by  Andersen  [Ref.  6]  and  Lee  [Ref.  7]. 
Finally  in  Section  D,  we  examine  120  absorbing  screens  of  uniform  heights  and  equal 
spacing  and  measure  the  relative  propagation  factor  at  the  final  range.  Again,  these 
screens  represent  a  row  of  buildings  or  houses  in  residential  areas  of  a  city.  The  result  of 
the  3D  PE  model  is  compared  with  the  numerical  integration  technique  result  proposed 
by  Bertoni  [Refs.  10]. 

A.  3D  PE  AND  THE  TWO-RAY  MODEL  RESULTS  COMPARISON 


In  this  section  we  stimulate  the  two-ray  model  and  the  3D  PE  model  with  the 
same  parameters  and  compare  their  results.  Figure  13  shows  the  two-ray  model  over  flat 
earth.  The  two-ray  model  propagation  factor,  F2ray,  [Ref.  2]  is  defined  as  follow  with 


2  ray 


\  +  Ye 


ik0AR 


(3.1) 


2H  H  sin  y/-Z  _i 

where  A R  = - ? ,  T  =  — — — ^ ,  y/ =  tan 


d 


sin  y/  +  Zs 


f  u  u  ^  _  18 cr(mS/m) 

,  +  ■ 


H,Hr 


V  d  j 


f(Mhz) 


and  Z,  =  7^— cos \y/)  /  £rc  -  -^=  \erc  \  » 


0. 
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The  H,  and  H,  are  the  heights  of  the  transmitter  and  receiver’s  antennas, 
respectively,  and  d  is  the  horizontal  distance  between  the  transmitter  and  receiver.  A R  is 
the  path  difference  between  direct  and  reflected  rays,  and  r  is  the  reflection  coefficient 
for  parallel  polarization  as  a  function  of  the  grazing  angle  \|/  and  the  complex  impedance 
Zs.  £rc  and  S  are  the  complex  dielectric  constant  and  conductivity  of  the  medium, 
respectively.  It  has  been  assumed  that  d  »  //,  and  Hr . 


Figure  13.  Two-ray  Model  Over  Flat  Earth. 

Figure  14  shows  the  vertical  cut  of  the  relative  propagation  factors,  F’s,  of  the 
two-ray  model  and  3D  PE  model.  The  two  models  are  stimulated  by  the  same  parameters 
as  listed  on  Figure  14.  The  results  were  taken  at  a  range  of  1000  m  from  the  transmitter  a 
over  flat  plane.  The  circle  with  the  solid  line  corresponds  to  the  two-ray  model  and  the 
plain  solid  line  is  the  3D  PE  result.  The  difference  between  the  two  results  is  within  1%. 
If  we  increase  the  height  of  the  receiving  antenna,  the  reflected  ray  will  tend  to  decrease 
and  the  direct  ray  will  dominate.  Therefore,  we  expect  the  two  F’s  to  approach  zero  dB 
for  large  heights. 
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Figure  14.  3D  PE  and  Two-ray  Comparison  (Vertical  Cut). 

B.  3D  PE  AND  FOUR-RAY  RESULTS  COMPARISON 

In  Section  A,  we  have  shown  that  the  result  of  the  3D  PE  agreed  with  the  two-ray 
results.  Now,  we  place  a  50  m  by  49.5  m  absorbing  screen  (a  single  knife-edge)  in  the 
3D  PE  model  located  125  m  from  the  transmitter  and  375  m  to  the  receiver  with  the 
frequency  of  operation  of  1  GHz.  The  loss  in  this  case  depends  on  the  height  of  the 
knife-edge  above  ground,  its  relative  location  from  the  transmitter/receiver,  the  ground 
constants,  and  the  frequency  of  operation.  Figure  15  shows  an  absorbing  knife-edge 
between  the  transmitter  and  the  receiver.  Table  2  provides  the  quantities  required  for  the 
evaluation  of  the  propagation  factor  Fk  e  for  the  2D-four-ray  model  [Ref.  2] .  The  3D  PE 

model  and  the  four-ray  model  are  simulated  for  the  same  assumed  parameters.  The  four 
paths  may  be  identified  from  the  transmitter  to  receiver  via  the  tip  of  the  knife-edge.  The 
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total  received  signal  is  the  sum  of  the  direct  wave,  if  any,  and  the  diffracted  waves 
received  via  each  one  of  the  four  paths.  The  total  diffracted  field  in  the  presence  of  the 
knife-edge  is  defined  by  equation  (3.2). 

(3-2) 

n=\ 


Figure  15.  An  Absorbing  Knife-Edge  Between  Transmitter  and  Receiver. 


Case  n 

Transmitter 

Location 

Receiver 

Location 

Free-Space  Path 
Length  rn 

Clearance 
Height  hn 

Reflection 
Coefficient,  rn 

1 

Tx 

Rx 

TxRx=p2+{H,-Hr)2 

//,  A +  ////, 

d 

1 

2 

Tx  ’ 

Rx 

TRx  =  sjd2  +  (Ht  +  Hr  )2 

H  +H,d2-Hrdl 
d 

Fa 

3 

Tx 

Rx’ 

TXRX  =  Txrx 

H  ha-ha 

k  d 

rB 

4 

Tx’ 

Rx’ 

TRX=TXRX 

ha+ha 

H‘+  d 

rArB 

Table  2.  Quantities  For  The  Evaluation  of  Fk.e{ After:  Ref.  2] 
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where  Eon  is  the  free-space  field  for  the  nth  situation,  Fn  is  the  normalized  field  for  that 
path,  and  H \  is  the  height  of  the  knife-edge.  Assuming  omni-directional  patterns  for  the 
transmitting  and  receiving  antennas  the  free-space  field  is  defined  as  follows 

„~jkorn 

Eon= - ,  «=1,...,4,  (3.3) 

r 


where  rn  is  the  free-space  path  length  for  the  nh  situation,  and  the  individual  knife-edge 
normalized  fields  are  given  by 


(1  +  J)F 


H  i 


r„, 


(3.4) 


where  F  is  the  Fresnel  integral,  rn  is  the  total  reflection  factor,  and  h„  is  the  clearance 


height  for  the  nh  path.  H /  is  the  radius  of  the  first  ellipsoid  in  the  plane  perpendicular  to 
the  line-of-sight  (LOS)  path  and  can  be  obtained  from  equation  (3.5): 


H  i  = 


\nXdxd2 
dx  +  d2 


(3.5) 


where  di  is  the  distance  between  the  transmitter  to  the  knife-edge,  dj  is  the  distance 
between  the  knife-edge  and  the  receiver,  and  X  is  the  operating  wavelength.  The 
propagation  factor  Fk  e  of  the  four-ray  model  can  then  be  obtained  from  equation  (3.6). 


k.e. 


(3.6) 


Additional  details  of  the  four-ray  model  are  provided  in  reference  2. 

Figure  16  shows  the  vertical  cut  of  the  relative  propagation  factors,  F’s,  of  the 
four-ray  model  and  3D  PE  model.  The  circle  with  the  solid  line  corresponds  to  the  four- 
ray  model  and  the  plain  solid  line  is  the  3D  PE  result.  The  results  were  taken  at  a  range 
of  500  m  from  the  transmitter  over  flat  plane.  The  screen  was  placed  125  m  from  the 
transmitter  and  375  m  to  the  receiver.  The  two  results  are  in  excellent  agreement.  As 
previously  mentioned,  the  relative  propagation  factors  F’s  approach  zero  dB  as  we 
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increase  the  receiving  antenna  height.  In  this  case,  it  is  already  close  to  zero  dB  at  a 
receiving  antenna  height  of  only  80  meters. 


Figure  16.  3D  PE  and  Four-ray  Comparison  (Vertical  Cut). 

C.  3D  PE  MODEL  AND  ANDERSEN’S  AND  LEE’S  RESULTS 
COMPARISON 

In  the  previous  section,  we  have  demonstrated  that  the  results  of  3D  PE  model 
corresponded  well  with  the  results  of  the  four-ray  model.  In  this  section,  we  consider 
nine  absorbing  screens  (multiple  knife-edges)  of  uniform  heights,  equal  spacing,  and 
variable  widths.  The  results  of  the  3D  PE  model  are  compared  with  the  results  presented 
by  Andersen  [Ref.  6]  and  Lee  [Ref.  7].  Andersen’s  results  were  obtained  from  the 
uniform  theory  of  diffraction  (UTD)  method.  Lee  presented  the  theoretical  results  in 
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reference  7.  Both  the  UTD  and  theoretical  results  assumed  that  the  screen  widths  are 
infinitely  long,  and  the  screen  heights  are  finite.  Full  details  on  the  UTD  method  are 
provided  by  references  6,  7,  and  14. 

We  also  examine  the  effects  of  three  finite  screen  widths  on  the  relative 
propagation  factors,  F’s.  The  three  finite  screen  widths  are  50  m,  25  m,  and  12.5  m.  In 
each  case,  nine  absorbing  screens  were  placed  between  the  transmitter  and  the  receiver 
with  each  screen  having  a  height  of  Hk  =10  m.  The  screens  are  placed  100  m  apart  from 
the  each  others.  The  first  screen  is  100  m  from  the  transmitting  antenna,  and  the  last 
screen  is  100  m  to  the  receiving  antenna.  The  distance  between  the  transmitter  and  the 
receiver  is  1000  m.  The  propagation  factors  are  determined  at  the  top  of  the  screens. 


The  theoretical  propagation  factor,  Fth,  in  dB  is  defined  by  equation  (3.7).  It  is  a 
function  of  edge  numbers,  where  the  edge  number  is  N  +1,  and  N  is  a  number  of  screens. 

f  1  ^ 


20  log 


10 


N  +  l 


(3.7) 


Again,  these  screens  represent  a  row  of  buildings  or  houses  in  residential  areas  of 
a  city.  In  these  types  of  environments,  the  base  station  antennas  of  cellular 
communication  systems  are  typically  located  above  or  near  to  the  rooftops  of  the 
surrounding  buildings.  In  these  cases,  the  propagation  takes  place  over  the  buildings, 
which  can  be  modeled  by  multiple  forward  diffraction  past  rows  of  buildings.  We  model 
each  row  of  buildings  as  an  absorbing  knife-edge,  and  via  the  3D  PE  numerical  technique 
we  determine  the  loss  associated  with  multiple  forward  diffraction  over  the  knife-edges. 
Furthermore,  depending  on  their  construction,  buildings  may  have  a  flat  roof,  a  peaked 
roof,  a  flat  roof  with  a  parapet,  or  a  myriad  of  other  roof  designs  [Ref.  13].  Figure  17 
illustrates  building  profiles  that  mobile  communication  systems  designers  can  be 
expected  to  encounter  for  typical  urban  environments  between  the  base  stations  and 
mobile  units.  Figure  17  (a)  shows  that  two  absorbing  screens  model  a  double  parapets 
roof  building.  Figure  17  (b)  illustrates  that  a  single  parapet  roof  building  is  modeled  by 
one  absorbing  screen  with  the  screen  placement  at  the  start  of  the  building.  And  Figure  17 
(c)  shows  that  a  peaked  roof  building  is  also  modeled  by  one  screen,  but  the  screen 
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placement  is  at  the  center  of  the  building  corresponding  to  the  peak  of  the  building. 
Reference  13  provides  more  details  on  building  profiles  and  screen  placements. 


Figure  17.  Typical  Building  Profiles  in  Urban  Areas,  with  Their  Equivalent  Screens 

Placements  [After:  Ref.  13]. 

Figure  18  illustrates  the  multiple  absorbing  knife-edges  with  equal  heights  and 
equal  spacing  for  this  test  problem.  The  assumed  parameters  used  in  the  3D  PE  model 
are  listed  on  Figure  19. 
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Figure  18.  Perfectly  Absorbing  Equal  Heights  and  Equal  Spacing  Multiple  Screens. 


Figure  19.  3D  PE,  UTD,  and  Theoretical  Results. 
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Figure  19  shows  the  results  of  the  3D  PE  model  with  the  results  presented  by 
Andersen  and  Lee.  The  solid  line  is  the  theoretical  result  for  2D.  The  circle  with  solid 
line  represents  the  results  presented  by  Andersen.  The  triangle  with  dashed  line  are  the 
results  of  50  m  screen  width,  the  diamond  with  solid  line  are  the  results  of  the  25  m 
screen  width,  and  the  star  with  dashed  line  are  the  results  of  thel2.5  m  screen  width.  The 
UTD,  the  50  m  and  25  m  screen  width  cases,  the  curve  slope  track  the  theoretical  closely 
up  to  the  seventh  screen.  After  the  seventh  screen,  they  are  only  0.5  dB  different  from 
the  theoretical.  For  the  case  of  12.5  m  screen  width,  the  results  do  not  quite  agree  with 
the  theoretical,  UTD,  50  m,  and  25  m  results,  but  the  difference  is  still  less  than  1  dB. 
These  differences  are  caused  by  the  fact  that  the  smaller  screen  width  allows  the  lateral 
waves  and  diffracted  waves  to  constructively  add  or  destructively  add.  These  lateral  and 
diffracted  waves  might  have  added  constructively  when  they  arrive  at  the  receiving 
antenna  up  to  the  seventh  screen;  then  added  destructively  after  the  seventh  screen,  which 
resulted  in  a  1  dB  different  from  the  theoretical  and  approximately  1.5  dB  different  from 
the  UTD  and  the  two  cases  of  the  3D  PE  model.  However,  we  expect  to  see  the  results  of 
the  3D  PE  model  approach  the  theoretical  results  if  we  consider  wider  screen  widths  that 
allow  less  lateral  waves  to  contribute  to  the  overall  results. 

All  presented  data  used  absorbing  screens  of  uniform  heights  and  equal  spacing, 
but  the  3D  PE  algorithm  is  capable  of  supporting  multiple  screens  of  variable  heights, 
widths,  and  spacing.  Based  on  the  comparison  results  in  sections  A,  B,  and  this  section, 
we  may  assume  that  the  3D  PE  model  is  capable  of  computing  the  propagation  factors 
over  flat  earth  with  multiple  absorbing  screens  of  non-uniform  heights,  unequal  spacing, 
and  variable  screen  widths.  Next  we  consider  120  screens  of  uniform  heights  and  equal 
spacing  over  flat  earth  and  perfectly  reflecting  ground. 

D.  3D  PE  MODEL  AND  120  UNIFORM  HEIGTH  AND  EQUAL  SPACING 

SCREENS 

In  this  section  we  consider  120  screens  of  uniform  heights  and  equal  spacing. 
Again,  the  propagation  takes  place  over  the  buildings.  As  previously  mentioned,  the  3D 
PE  model  is  capable  of  supporting  multiple  absorbing  screens  of  variable  heights,  widths, 
and  spacing.  But  for  the  purpose  of  this  study,  we  will  only  consider  absorbing  screens 
of  uniform  heights  and  equal  spacing.  We  compare  the  3D  PE  model  results  with  the 
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results  presented  by  Bertoni  [Ref.  10].  Figure  20  shows  setup  geometry  of  this  test  case. 
We  consider  a  plane  wave  incident  at  an  angle  a  =  1°  at  an  operating  frequency  of  900 
MHz.  The  distance  between  screens  is  50  m,  which  means  that  the  distance  between  the 
transmitter  and  the  receiver  is  6050  m  apart.  Based  on  the  incident  angle  a  of  1°  and  the 
separation  distance  between  the  transmitter  and  the  last  screen,  we  determine  that  the 
vertical  clearance  height,  Hch,  from  the  top  of  the  120th  screen  to  the  tip  of  the 
transmitting  antenna  is  105  m.  If  the  screen  height  of  20  m  is  chosen,  then  the 
transmitting  antenna  height  should  be  125  m.  We  assume  the  screens  width  to  be  50  m. 


Figure  20.  120  Equal  Height  and  Equal  Spacing  Screens. 

Figure  21  shows  the  height  variation  of  the  relative  field  strength,  F,  computed  by 
3D  PE  model,  incident  on  the  row  of  120  screens  of  20  m  heights  and  50  m  spacing  for 
the  parameters  assumed  the  previous  paragraph.  The  receiving  antenna  height  is 
measured  in  wavelength  (m)  while  the  magnitude  of  the  field  strength  is  in  linear  scale. 
Figure  19  also  shows  the  simple  diffraction  in  the  shadow  region  of  the  receiving  antenna 
H,  <  0.  We  use  the  rooftop  as  the  reference  point  for  H,  =  0.  Above  the  rooftops,  the 
field  variation  is  similar  to  that  of  a  standing  wave  resulting  from  the  summation  of  the 
direct  waves  and  reflected  waves  from  the  plane  of  the  rooftops.  The  reflected  waves  are 
in  fact  the  sum  of  the  waves  diffracted  from  the  rooftops,  whose  phase  variations  cause 
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them  to  add  constructively  in  the  specular  direction  and  other  grating  lobe  directions. 
Because  the  diffraction  coefficients  decrease  with  angle,  the  reflected  field  is  greatest  in 
the  specular  direction. 


Ny  =  1 024  Nz  =  1 024  Dx  =  25  meters  Ht  =  1 25  meters 


Receiving  Antenna  in  Wavelegth 


Figure  21.  120  3D  PE  Results  For  Equal  Height  and  Equal  Spacing. 

Figure  22  shows  the  height  dependence  of  the  field  strength,  F,  of  the  3D  PE 

model  with  the  result  from  the  numerical  integration  method  presented  by  Bertoni  and  the 

result  from  the  UTD  [Ref.  10].  The  circle  with  solid  line  is  the  sampling  result  of  the  3D 

PE  model  shown  in  Figure  21.  The  asterisk  with  solid  line  is  the  result  of  the  numerical 

integration  method,  and  the  triangle  with  solid  line  is  the  result  of  the  UTD.  In  the 

shadow  region,  the  results  of  the  3D  PE  model  and  the  numerical  integration  method 

have  good  agreement,  but  the  result  of  the  UTD  is  slight  higher,  which  implies  that  the 

UTD  overestimates  the  field  amplitude  [Ref.  10].  Above  the  shadow  region,  the 
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IV.  HILLY  TERRAINS  RESULTS 


In  this  chapter,  the  propagation  takes  place  over  hilly  terrains,  and  we  study  the 
field  strengths.  A  single  round  is  considered  in  Section  A.  The  two  hills  of  sinusoidal 
shape  are  considered  in  Section  B.  In  both  cases,  we  place  multiple  absorbing  screens  of 
uniform  heights,  equal  spacing,  and  variable  widths  between  the  transmitter  and  the 
receiver.  We  measure  the  field  strengths,  instead  of  propagation  factors,  at  the  rooftops. 
The  results  of  the  3D  PE  model  are  compared  with  the  results  presented  by  Piazzi  and 
Bertoni  [Refs.  8  -  10]. 

A.  A  ROUND  HILL  AND  3D  PE  MODEL 

In  the  previous  test  cases,  we  have  examined  the  propagation  factors  over  flat 
earth  with  multiple  screens  of  uniform  heights,  equal  spacing,  and  variable  widths.  Now, 
we  consider  a  rounded  hill  with  multiple  absorbing  screens  of  uniform  heights,  equal 
spacing,  and  variable  widths.  We  examine  the  diffracted  field  strengths  at  the  rooftops, 
instead  of  the  propagation  factors  as  in  the  previous  cases.  We  also  consider  the  effects 
of  finite  screen  widths  on  the  field  strengths.  We  compare  simulation  results  of  the  3D 
PE  model  with  results  presented  by  Piazzi  [Ref.  8].  Figure  23  illustrates  setup  geometry 
for  this  test  case.  These  screens  represent  a  row  of  buildings  or  houses  in  residential 
areas  of  a  well  built-up  city.  The  buildings  are  7  m  high  and  the  row  separation  Ds  is  50 
m,  and  the  screen  widths  are  Wk  =  100  m,  50  m,  and  25  m.  The  cylindrical  hill  has  a 
radius  R/,  of  10  km  [Ref.  8].  The  maximum  height  of  the  hill  is  50  m.  The  base  of  the 
hill  is  at  a  distance  of  X/,  =1000  m  from  the  peak.  For  the  simulation,  the  transmitting 
antenna  is  assumed  to  be  located  at  the  point  x  =  -1000  m  from  the  hilltop  and  has  a 
height  of  H,  =  57  m,  which  places  it  at  the  same  height  as  the  highest  rooftops  [Ref.  10]. 
The  current  source  is  excited  at  an  operating  frequency  of  900  MHz.  The  rounded  hill  is 
defined  by  equation  (4.1). 

Z  =  (4.1) 

where  Ri,  is  the  radius  of  the  hill,  Hi  is  the  height  of  the  buildings,  and  x  is  the  range  in 
meter.  The  total  simulation  range  is  6000  m. 


41 


Y 


Figure  23.  A  Round  Hill  with  Multiple  Screens  of  Uniform  Heights  and  Equal 

Distances. 

Figure  24  shows  the  profile  of  multiple  absorbing  screens  representing  Figure  23 
used  to  carry  out  the  simulation  in  the  3D  PE  model.  Furthermore,  if  we  can  model  an 
urban  area  with  multiple  absorbing  screens  of  variable  heights,  widths,  and  spacing,  the 
3D  PE  model  can  evaluate  the  propagation  factors  as  well  as  field  strengths  at  any  point 
in  the  simulation. 

Figure  25  shows  the  results  of  3D  PE  model,  the  free-space  result,  and  Piazzi’s 
result.  There  are  three  3D  PE  model  results.  The  solid  line  is  the  result  of  the  3D  PE 
model  with  100  m  screen  width.  The  dashed  line  is  the  result  of  the  50  m  screen  width. 
The  dashed  and  dot  line  is  the  result  of  the  25  m  screen  width.  The  triangle  with  solid 
line  is  the  free-space  result,  and  the  circles  are  Piazzi’s  results.  Up  to  1000  m  away  from 
the  base  station,  the  field  strength  behaves  like  free-space.  Then  at  the  hilltop  the  field 
behaves  like  a  single  knife-edge  case,  where  the  magnitude  of  the  field  strength  is  6  dB 
below  the  free-space.  Between  the  base  station  and  the  peak  of  the  hill,  the  field  strength 
does  not  behave  like  a  free-space  field.  Perhaps  this  is  one  of  the  limitations  for  the  3D 
PE  model,  but  this  is  minor  issue.  After  the  hilltop  and  behind  the  hill  up  to  x  =  1000  m, 
the  field  strength  decays  rapidly  and  linearly  with  distance  at  a  rate  of  approximately 
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0.055  dB  per  meter.  This  may  be  understood  in  terms  of  the  creeping  wave 
representation  of  the  fields  traveling  over  the  houses.  At  the  point  of  x  =  1000  m  from 
the  hilltop,  the  field  strength  reaches  its  lowest  point,  but  begins  to  increase  rapidly  to 
maximum  and  then  decreases  slowly  with  distance.  We  expect  the  field  to  continue 
decaying  as  it  propagates  away  from  the  base  of  the  hill,  but  that  is  not  the  case  because 
the  diffracted  waves  from  the  rooftops  on  the  hillside.  We  also  expect  the  field  to 
eventually  level  off  and  begin  to  decay  and  behave  like  multiple  knife-edges  propagation 
case.  As  shown  in  Figure  25,  the  field  levels  off  around  x  =  3000  m  away  from  the  base 
of  the  hill.  If  we  continue  to  measure  the  field  strength  for  both  the  3D  PE  model  and  the 
free-space,  the  field  strength  of  the  3D  PE  model  will  eventually  be  6  dB  below  the  free- 
space  value  as  in  the  a  single  knife-edge  case. 


Figure  24.  Side  View  of  a  Single  Hill  with  Rows  of  Uniform  Heights  and  Spacing 


Buildings. 

43 


Figure  25  also  shows  that  we  have  excellent  agreement  between  Piazzi’s  result 
and  the  result  of  the  3D  PE  model  especially  with  screen  width  Wk  =100  m.  However,  it 
is  not  the  case  for  the  50  m  and  25  m  screen  widths.  Furthermore,  for  both  the  50  m  and 
25  m  screen  widths,  there  are  oscillations  as  the  fields  move  away  from  the  base  of  the 
hill.  There  appears  to  be  more  oscillation  in  the  25  m  screen  width  case  than  for  the  50  m 
case.  Its  magnitude  is  also  higher  than  the  other  cases,  but  it  begins  to  approach  other 
cases  as  it  moves  away  from  the  base  of  the  hill  on  backside.  The  oscillations  and  higher 
magnitudes  are  caused  by  lateral  waves  because  of  the  finite  screen  widths.  Therefore,  if 
we  increase  the  screen  width,  then  we  expect  the  result  3D  PE  model  to  agree  with 
Piazzi’s  result  because  Piazzi  assumed  the  screen  width  is  infinitely  long.  This  is 
equivalent  to  the  2D  PE  model  which  did  not  take  into  account  the  lateral  waves.  This  is 
demonstrated  the  100  m  screen  width  case,  which  prevented  the  lateral  waves  from 
significantly  contributing  to  the  overall  result.  We  have  shown  that  the  3D  PE  model  is 
capable  of  supporting  the  variable  screen  widths.  The  last  test  problem  we  will  present  in 
this  thesis  is  the  two  sinusoidal  hills. 


Screen  Position  (m) 


Figure  25.  A  Single  Hill  with  Rows  of  Buildings  of  Uniform  Heights  and  Spacing 

Field  Strength. 
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B.  3D  PE  MODEL  AND  TWO  HILLS  OL  SINUSOIDAL  SHAPE 


In  the  previous  section,  we  have  considered  a  single  rounded  hill  and 
demonstrated  that  the  result  of  the  3D  PE  model  agreed  with  Piazzi’s  result.  We  have 
also  shown  that  the  3D  PE  model  can  support  variable  screen  widths.  Now  we  consider 
two  sinusoidal  hills,  and  this  is  the  last  test  problem  we  will  consider  for  this  thesis.  It 
has  similar  setup  geometry  as  shown  in  Figure  23,  except  it  has  two  hills  of  sinusoidal 
shape.  Figure  26  illustrates  setup  geometry  for  this  test  case. 


Figure  26.  Two  Hills  of  Sinusoidal  Shape. 


Again,  these  screens  represent  a  row  of  buildings  or  houses  in  residential  areas  of 
a  well  built-up  city.  The  buildings  are  7  m  high  and  the  row  separation  Ds  is  50  m  with 
building  widths  W*  =  100  m,  50  m,  and  25  m.  The  maximum  heights  of  the  two  hills  are 
50  m.  The  two  sinusoidal  shape  hills  have  period  of  3000  m  [Ref.  9],  which  means  the 
distance  from  peak  to  peak  is  3000  m.  The  distance  A/,  from  the  first  peak  to  its  base, 
measured  toward  the  second  peak,  is  1500  m.  For  the  simulation,  the  transmitting 
antenna  is  assumed  to  be  located  at  the  point  x  =  -1500  m  from  the  first  hilltop  and  has  a 
height  of  Ht  =  57  m,  which  again  places  it  at  the  same  height  as  the  highest  rooftops  [Ref. 
10].  The  two  hills  of  sinusoidal  shape  are  defined  by 


Z(x)  =  Hk  +25  +  25cos 


^  2nx 
v  3000 , 


(4.2) 
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where  H \  is  the  height  of  the  buildings  and  x  is  the  range  in  meters. 

Figure  27  shows  the  profile  of  multiple  absorbing  screens  representing  Figure  24 
used  to  carry  out  the  simulation  in  the  3D  PE  model.  The  operating  frequency  is  900 
MHz  and  same  as  the  previous  case. 


Screen  Placement  (m) 


Figure  27.  Side  View  of  two  Sinusoidal  Hills  plus  Buildings  Profile. 

Figure  28  shows  the  composite  results  of  the  3D  PE  model,  free-space  and  Piazzi. 
The  solid  line  is  the  result  of  the  3D  PE  model  with  100  m  screen  width.  The  dashed  line 
is  the  result  of  the  50  m  screen  width.  The  dashed  and  dot  line  is  the  result  of  the  25  m 
screen  width.  The  triangle  with  solid  line  is  the  free-space  result.  And  the  circles  are 
Piazzi’s  results.  Up  to  1500  m  away  from  the  base  station,  the  field  strength  behaves  like 
free-space,  which  corresponds  to  the  peak  of  the  first  hill.  Then  at  the  first  hilltop  the 
field  behaves  like  a  single  knife-edge  problem,  same  as  the  previous  case,  where  the 
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magnitude  of  the  field  strength  is  6  dB  below  the  free-space.  As  in  the  one  hill  problem, 
the  field  does  not  behave  like  free-space  field  when  it  is  close  to  the  base  station.  The 
field  strength  minima  occur  at  x  =  1000  and  4000  m  which  do  not  correspond  to  the  bases 
of  the  hills  as  it  did  in  a  single  round  hill  case.  When  we  decrease  the  widths  of  the 
screen,  we  see  more  oscillations  as  demonstrated  in  the  one  hill  case.  The  oscillation 
appears  to  be  more  severe  in  the  25  m  screen  width  case.  These  oscillations  are  caused 
by  the  lateral  waves. 


Figure  28.  Two  Sinusoidal  Hills  with  Rows  of  Buildings  of  Uniform  Heights  and 

Spacing  Field  Strength. 

We  have  thus  demonstrated  that  the  3D  PE  model  can  support  both  flat  and  hilly 
terrains  with  multiple  absorbing  screens  of  uniform  heights,  equal  spacing,  and  variable 
screen  widths. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

A  3D  model  based  on  the  parabolic  equation  (PE)  for  predicting  propagation  path 
loss  has  been  developed  in  this  thesis.  Two  types  of  terrains  were  considered:  the  flat 
earth  and  the  hilly  terrain.  For  the  flat  earth  case,  four  test  problems  were  examined.  The 
results  from  these  cases  were  compared  with  results  available  from  the  two-ray  model, 
the  four-ray  model,  the  uniform  theory  of  diffraction  (UTD),  Lee’s  formulation,  and  the 
numerical  integration  technique  as  proposed  by  Bertoni.  For  the  hilly  terrain  case,  two 
test  problems  were  considered.  The  results  from  this  case  were  compared  with  the  results 
presented  by  Piazzi  and  Bertoni.  For  all  these  test  problems,  ground  was  assumed  to  be 
non-perfectly  reflecting,  and  buildings  as  perfectly  absorbing  because  all  available  results 
that  were  used  in  the  comparisons  made  these  assumptions. 

For  the  flat  earth  case,  the  first  test  problem  was  to  simulate  and  determine  the 
propagation  factor  F  at  the  final  range  over  flat  earth  without  any  obstacles  between  the 
transmitter  and  the  receiver.  The  transmitting  antenna  had  a  height  H,  =  30  m,  and  the 
operating  frequency  was  1  GHz.  The  receiver  located  1000  m  from  the  transmitter.  The 
results  of  the  3D  PE  model  compared  excellently  with  those  of  the  two-ray  model. 

The  second  test  problem  was  to  simulate  and  determine  F  with  a  single  absorbing 
screen  placed  between  the  transmitter  and  the  receiver  that  had  a  height  =  50  m,  and  a 
width  Wk  =  49.5  m.  The  screen  representing  a  building  was  located  125  m  from  the 
transmitter  and  375  m  to  the  receiver.  The  transmitting  antenna  height  H,  =  60  m,  and  the 
operating  frequency  was  1  GHz.  The  results  of  the  3D  PE  model  in  this  case  were 
compared  with  the  results  of  the  four-ray  model.  Once  again,  the  results  had  excellent 
agreement.  F  approached  0  dB  as  the  receiving  antenna  height  increased.  In  this  case, 
both  the  3D  PE  model  and  the  four-ray  results  approached  0  dB  at  the  receiving  antenna 
height  of  80  m. 

In  the  third  test  problem,  nine  absorbing  screens  of  uniform  heights  and  equal 
spacing  were  placed  between  the  transmitter  and  the  receiver.  This  test  case  was 
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equivalent  to  a  radiowave  propagating  over  a  row  of  houses  in  an  urban  area.  The 
screens  had  a  height  Hi  =  10  m  and  the  widths,  W/(.  The  three  different  screen  widths  Wi_ 
were  12.5  m,  25  m,  50  m.  Each  screen  had  a  height  of  10  m  and  a  spacing  distance  of  50 
m.  Their  effects  on  the  propagation  factor  F  were  studied.  The  propagation  factor  in  this 
case  was  determined  at  the  rooftops.  The  results  of  the  3D  PE  model  were  compared 
with  the  UTD  results  presented  by  Andersen  and  the  theoretical  model  results  proposed 
by  Lee.  The  3D  PE  model  results  with  the  screen  widths  of  25  m  and  50  m  had  excellent 
agreement  with  the  UTD,  and  both  the  UTD  and  the  3D  PE  model  results  followed  the 
theoretical  results  only  up  to  the  seventh  screen.  Then  up  to  the  ninth  screen  they  were 
0.5  dB  higher  than  the  theoretical  predictions.  For  the  12.5  m  case,  the  results  were 
slightly  higher  than  all  cases  up  to  the  seventh  screen;  then  they  were  lower  than  all  cases 
up  to  the  ninth  screen.  These  differences  were  caused  by  the  lateral  waves  due  to  the 
finite  screen  width. 

The  last  test  problem  for  the  flat  earth  case  involved  120  absorbing  screens  of 
uniform  heights  and  equal  spacing  between  the  transmitter  and  the  receiver.  These 
screens  represented  row  of  buildings  or  houses.  The  transmitting  antenna  had  a  height  of 
H,  =  125  m,  and  the  operating  frequency  was  900  MHz.  The  screens  had  the  height,  H/c  = 
20  m,  and  the  width,  W*  =  50  m.  They  were  spaced  50  m  apart.  The  propagation  factor 
was  determined  at  the  final  range.  The  results  were  compared  with  the  results  of  the 
numerical  integration  technique  presented  by  Bertoni  and  the  UTD.  In  the  shadow 
region,  the  results  of  the  3D  PE  and  the  numerical  integration  method  agreed,  but  the 
result  of  the  UTD  was  slightly  higher,  which  implied  that  the  UTD  overestimated  the 
field  amplitude.  Above  the  shadow  region,  the  numerical  integration  method  and  the 
UTD  had  excellent  agreement,  and  although  the  results  of  the  3D  PE  model  were  0.25  dB 
lower  than  those  two  models,  nonetheless  there  was  reasonably  good  agreement  with 
them. 

For  the  hilly  terrain  case,  the  first  test  problem  was  a  single  rounded  hill  with 

multiple  absorbing  screens  of  uniform  heights  and  equal  spacing.  The  screens  had  the 

heights  Hk  =  7  m.  The  field  strengths  were  determined  on  rooftops,  and  were  compared 

with  the  numerical  integration  technique  results  presented  by  Piazzi  and  Bertoni.  The 

second  test  problem,  the  propagation  took  placed  over  two  hills  of  sinusoidal  shape  with 
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multiple  absorbing  screens  of  uniform  heights  and  equal  spacing.  Again,  the  results  were 
compared  with  the  results  presented  by  Piazzi  and  Bertoni.  The  screens  in  this  case  also 
had  the  heights  Hk  =  7  m  In  both  cases,  three  variable  screen  widths  Wk  were  considered: 
25  m,  50  m  and  100  m.  Their  effects  on  the  field  strengths  were  studied.  For  both  cases, 
the  transmitting  antenna  height  was  H,  =  57  m,  and  the  operating  frequency  was  900 
MHz. 

For  a  single  round  hill,  the  result  of  the  3D  PE  model  with  screen  width  of  100  m 
compared  excellently  with  the  Piazzi’ s  result.  Before  the  hilltop,  the  field  strength 
behaved  like  free-space.  Then  at  the  hilltop  the  field  behaved  as  in  a  single  knife-edge 
case,  where  the  magnitude  of  the  field  strength  was  6  dB  below  the  free-space.  On  the 
backside  of  the  hill,  the  field  strength  decayed  rapidly  and  linearly  with  distance,  but 
began  to  increase  rapidly  to  maximum  and  then  decreased  slowly  with  distance.  For  the 
cases  of  25  m  and  50  m  screen  width,  the  field  strengths  agreed  with  the  100  m  case  and 
Piazzi’ s  result  until  after  the  base  of  the  hill  on  the  backside  after  which  it  began  to 
increase  rapidly  with  oscillations.  These  oscillations  were  caused  by  the  laterally 
diffracted  waves.  The  oscillations  were  more  prominent  in  the  25  m  case  than  the  50  m 
case. 

For  two  sinusoidal  hills,  the  field  strengths  behaved  like  free-space  up  to  the 
first  hilltop  similar  to  the  previous  case.  Then  at  the  first  hilltop  the  fields  behaved  like  a 
single  knife-edge  problem,  where  the  magnitude  of  the  field  strength  is  6  dB  below  the 
free-space.  The  field  strength  minima  occurred  at  x  =  1000  m  and  4000  m  which  did  not 
correspond  to  the  bases  of  the  hills.  When  the  widths  of  the  screen  were  decreased,  the 
oscillations  showed  up  similarly  to  the  single  rounded  hill  case.  The  oscillations 
appeared  to  be  more  severe  in  the  screen  width  of  25  m. 

Finally,  the  3D  PE  model  for  predicting  propagation  path  loss  in  urban  areas  on 
flat  and  hilly  terrains  was  developed.  Six  different  test  problems  were  considered.  The 
results  were  compared  with  the  results  available  in  open  literature  and  showed  excellent 
agreement.  We  have  also  demonstrated  that  the  3D  PE  model  can  support  both  flat  and 
hilly  terrains  with  multiple  absorbing  screens  of  uniform  heights,  equal  spacing  and 
variable  screen  widths. 
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B.  RECOMMENDATIONS 


Since  realistic  environments  will  have  non-perfectly  reflecting  obstacles,  a  more 
complete  3D  PE  formulation  for  non-perfectly  reflecting  obstacles  should  be  studied  as 
well  as  the  effects  of  the  antenna  depolarization.  In  addition,  it  is  recommended  that  the 
algorithm  be  implemented  in  C  or  Fortran  for  faster  computation,  along  with  the 
automating  terrain  and  building  inputs  for  the  model.  It  may  be  desirable  to  add  the 
capability  to  utilize  the  Digital  Terrain  Elevation  Data  (DTED)  and  Digital  Feature 
Attribute  Data  (DFAD)  available  at  the  National  Imagery  and  Mapping  Agency  (NIMA). 
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APPENDIX  A:  3D  PE  ALGORITHM  AND  2RAY  MODEL 


Below  is  the  MatLab  program  listing  of  the  3D  parabolic  equation  model  and  the 
2Ray  model.  This  program  will  compute  the  propagation  factors  for  the  3D  PE  with  the 
2Ray  model  over  perfectly  flat  earth  and  constant  admittance  plane.  The  operating 
frequency  is  1  GHz. 


o 

%  Declaring  Contants 

o 


hold  off; 
tic 
c=3e8 ; 
f=le9; 

Lamda=c/f ; 

ko= (2*pi) /Lamda; 

S=2  5 ; 

Er=5  ; 


%  Speed  of  light 
%  Operating  frequency 
%  Wavelength 
%  Wavenumber 

%  Ground  conductivity  in  mS/m 
%  Relative  dielectric  constant 


Dz=Lamda; 

Dy=Lamda; 

sigmaz=Lamda; 


Delta  Z 
Delta  Y 
Current 


in  meter 
in  meter 

source  standard  deviation 


9?iri?iriti?i?iririri?iri?iri?iritiri?iri?'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%  Input  Parameters 

o 


Ht=30;  %  Transmitting  antenna  height 


D=1000 ; 

o, 

o 

Range  from 

the 

transmitter 

Dx=100 ; 

o. 

o 

Incremental 

range 

(Dx)  in  meter 

Ny=1024 ; 

o, 

o 

Sample  size 

in 

the 

y-direction 

Nz=l 02  4 ; 

o. 

o 

Sample  size 

in 

the 

z-direction 

S^iritiritiritiririririririritiririritirir'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%  Define  range  in  y-direction  and  z-direction 

o 


ymax=Ny*Dy/ 2 ; 
yl=0 : Dy : ymax; 
y=[yl  zeros ( 1 , Ny/2-1 )] ; 


%  maximum  range  in  y-direction 

%  First  half  of  y  range  incrementally  increase 
%  y  array 


zmax  =  Nz*Dz/2; 

zl=0 : Dz : zmax; 

z— [ z 1  zeros (l,Nz/2-l) ] ; 


%  Find  maximum  z  range  (Zmax) 

%  First  half  of  the  range  of  Zmax+1 
%  Construct  a  full  z  array 
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o 

%  Define  Wavenumbers  in  Spatial  Domain 

o 


kymax=pi/Dy; 

Dky=2  *kymax/Ny ; 
kyl=-kymax : Dky : kymax; 
ky=[kyl(:, (Ny/2 ) +1 : Ny+1 )  kyl 


%  Maximum  wavenumber  in  y-direction 
%  detla  ky 
%  Range  of  ky 
, 2 :Ny/2) ] ; 


kzmax=pi/Dz ; 

Dkz=2  *kzmax/Nz ; 
kz l=-kzmax : Dkz : kzmax ; 
kz=[kzl(:, (Nz/2) +1 :Nz+l)  kzl 


%  Maximum  wavenumber  in  z-direction 
%  detla  kz 
%  Range  of  kz 
,2:Nz/2) ] ; 


o 

%  Compute  Wavenumbers  in  Frequency  Domain,  kx 

o 


Ky=meshgrid (ky,  1 

: Nz )  ; 

o, 

o 

Create  a  NyxNz  matrix 

of 

ky  row-repeat 

Kz=meshgrid (kz ,  1 
repeat 

:Ny) . ' ; 

o, 

o 

Create  a  NyxNz  matrix 

of 

kz  column- 

kx=sqrt (koA2-Ky. 

A2-Kz. A2) ; 

o, 

o 

Compute  theNy  x  Nz  kx 

mat 

rix 

clear  Ky; 

o, 

o 

Clear  matrices  Ky  and 

Kz 

from  memory 

clear  Kz; 

o, 

o 

To  make  the  algorithm 

run 

faster 

o 

%  Compute  reflection  coefficient 

o 


Erc=Er+i*18*S/ (f/le6) ; 
Zs=l/sqrt (Ere) ; 

Gamma= (kz-ko*Zs) ./ (kz+ko*Zs) ; 
%gamma=-ones (Nz,Ny) ; 
gamma=meshgrid (Gamma, 1 :Ny) . ' ; 


%  Complex  dielectric 
%  Impedance 

%  Complex  reflection  coefficient 
%  Setting  gamma  =  -1 

%  Create  reflection  coefficient  matrix 
%  by  column  repeating 


o 

%  Define  the  Hanning  Window 

o 


Hya= [ ] ;  % 

for  t=0:Ny/2;  % 

if  (t>=0  &  t<=3*Ny/8)  % 

hy=l ; 

elseif  (t>=3*Ny/8  &  t<=Ny/2)% 
hy= (sin (4*pi*t/Ny) ) A2; 

end 

Hya= [Hya  hy] ;  % 

end  % 

Yy=fliplr (Hya ( : , 2 :Ny/2) )  ;  % 

second  element  to 

o, 

o 

HY= [Hya  Yy] ;  % 


An  empty  vector 
Define  Ny/2  +  1  points 
Define  first  3Ny/8  +1  points 

Define  next  Ny/8  points 

Construct  the  first  half  of  Hanning 
window  of  l+Ny/2  elements 

Flip  left  to  right  of  the  Hanning 
%  windows  from  the 

Ny/2  element  for  total  of  [lx(Ny/2-l)] 
Hanning  window  in  y-direction 
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Hzb=  [  ] ;  % 

for  t=0:Nz/2;  % 

if  ( t >— 0  &  t<=3*Nz/8)  % 

hz  =  l ; 

elseif  (t>=3*Nz/8  &  t<=Nz/2)% 
hz=(sin(4*pi*t/Nz)  )  A2; 

end 

Hzb= [Hzb  hz] ;  % 

end  % 


An  empty  vector 
Define  Nz/2  +  1  points 
Define  first  3Nz/8  +1  points 

Define  next  Nz/8  points 

Construct  the  first  half  of  Hanning 
window  of  l+Nz/2  elements 


Yz=fliplr(Hzb(:,2:Nz/2) )  ; 


HZ= [Hzb  Yz] ' ; 


%  Flip  left  to  right  of  the  Hanning 
%  windows  from  the  second  element  to 
%  Nz/2  element  for  total  of  [lx(Nz/2-l)] 
%  Hanning  window  in  z-direction 


Hmy=meshgrid (HY, 1 : Nz )  ; 
Hmz=meshgrid (HZ , 1 : Ny )  '  ; 


%  Row  repeat 
%  Column  repeat 


H3D= (Hmy . *Hmz ) ; 


%  3D  Hanning  window 


o 

%  Gaussian  Current  Source 

o 


f= (1/ (sqrt (2*pi) *sigmaz) ) *exp (- (z-Ht) . A2/ (2*sigmazA2) ) ; 
f_tilda=Dz*f ft (f ) ;  %  Fourier  transform  of  f(z) 

o 

%  Initial  H-field  at  x=0 

o 

g_tilda=exp (-kz . A2*sigmazA2/2) ;  %  initial  g  tilda 

hyeO_tilda=g_tilda . *cos (kz*Ht) ;  %  Initial  hye_tilda (0+, ky, kz)  field 

hyoO_tilda=-i*g_tilda . *sin (kz*Ht) ;  %  Initial  hyo_tilda (0+, ky, kz)  field 

%  Initial  even  H-field,  column  repeat 
HyeO_tilda=meshgrid (hyeO_tilda, 1 :Ny) . 

%  Initial  odd  H-field,  column  repeat 
HyoO_tilda=meshgrid (hyoO_tilda, 1 :Ny) . 

%  Include  the  reflection  coefficient 
Hy0_tilda=0 . 5* (1+gamma) . *Hye0_tilda+0 . 5* (1-gamma) . *HyoO_tilda; 

Hy_tilda=HyO_tilda . *H3D; %  Hy_tilda (0+, ky, kz)  at  x=0 

Exp2=exp (i*kx*Dx) ;  %  Marching  range 
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o 


%  3D  Parabolic  Equation  Basic  Algorithm 

S^i?i?iri?i?i(i?i?i?i?iri?iri(i?i?i?itiri?ici(*i?i?iti?i?'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 


for  x=0:Dx:D-Dx; 

Hy_tilda_Dx=Hy_tilda . *Exp2 ; 


%  Hy_tilda (Dx,  ky ,  kz ) 


Hy=Dkz*Dky* (Ny*Nz*ifft2 (Hy_tilda_Dx) ) / (2*pi) A2; 

%  ifft2  wrt  to  ky  and  kz 


Hy_H=Hy . *H3D; 

Hy z 1 =Hy_H ( 1 : N z / 2  + 1 ,  :) ; 

HyzO=flipud (Hy_H (2 :Nz/2, : ) ) ; 

Hye=[Hyzl;  HyzO]; 

Hyo=[Hyzl;  -HyzO]; 

Hye_tilda= (Dz*Dy) *fft2 (Hye) ; 
Hye_tilda_H=Hye_tilda . *H3D ; 

Hyo_tilda= (Dz*Dy) *fft2 (Hyo) ; 
Hyo_tilda_H=Hyo_tilda . *H3D ; 


%  Apply  the  Hamming  window 
%  H-field  for  z  >  0 
%  H-field  for  z  <  0 
%  Even  H-field 
%  Odd  H-field 

%  Take  the  FFT  of  the  even  H~ 

%  Apply  the  Hanning  window  in  freq 

%  Take  the  FFT  of  the  odd  H~ 

%  Apply  the  Hanning  window  in  freq 


%  Multiply  by  the  reflection  coefficient 
Hye_tilda_g=0 . 5* ( 1+gamma) . *Hye_tilda_H; 

%  Multiply  by  the  reflection  coefficient 
Hyo_tilda_g=0 . 5* ( 1-gamma) . *Hyo_tilda_H; 

Hy_tildal=Hye_tilda_g  +  Hyo_tilda_g; %  Hy_tilda (x, ky , kz ) 
Hy_tilda=Hy_tildal . *H3D;  %  Apply  the  Hanning  window 


end 

o 

%  Compute  Propagation  Factor  at  Final  Range 

o 

J=2;  %  Pick  a  column  of  H-field  at  final  range 

rho=sqrt (DA2+ ( J) A2 )  ;  %  Distance  from  transmitter  to  receiver 

R=sqrt ( ( z ( : , 1 : Nz /2+1 ) -Ht ) . A2  +  rhoA2); 

%  Distance  from  transmitter  to  receiver 
theta=asin ( rho . /R) ;  %  Angle  Theta  in  spherical  coordinate 

phi=asin ( J/rho) ;  %  Angle  Phi  in  spherical  coordinate 

g_tilda_f s=exp (- ( (cos (theta) *ko) . A2*sigmazA2) /2) ; 

%  Free-space  g_tilda 

Hyfs=-i*ko*sin (theta) *cos (phi) . *g_tilda_fs . *exp (i*ko*R) ./ (4*pi*R) ; 

%  Free-space  H-field 

A=Hy_H (1 :Nz/2+l, 2) ;  %  z  >  0  of  a  column  of  the  H-field 

F=20*logl0 (abs (A. /Hyfs ' ) ) ;  %  Propagation  factor 

save  H:\THESIS_2001\MATLAB_CODES\PERF2Ray  F  -ASCII 
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o 


%  2Ray-Model 

S^*i?iri?iciticiti?i?iri?i?iti?itirifiri?iri(*i?i?i?'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 


CN 

II 

o, 

o 

rho=sqrt (DA2+ (J) A2) ; 

o, 

o 

receiver 

Hr=z (1, 1 :Nz/2+l) ; 

o, 

o 

DR=2*Ht*Hr/ rho; 

O. 

O 

Q. 

psi=atan ( (Ht+Hr) /rho) ; 

O 

O, 

O 

A  column  final  H-field 

Distance  from  transmitter  to  the 

The  receiving  antenna  height 
Path  difference  between  direct  and 
ground  reflected  waves 
Grazing  angle 


Z2ray=sqrt (Erc-cos (psi) . A2) /Ere; %  Normalized  surface  impedance  Z(psi) 
gamma2ray= ( (sin (psi) -Z2ray) . / (sin (psi) +Z2ray) ) ; %  Reflection  coefficient 
Tworay=abs (l+gamma2ray . *exp (i*ko*DR) ) ; %  2ray  model  propagation  factor 
F2ray=20*logl0 (Tworay) ;  %  2ray  model  propagation  factor  in  dB 

save  H:\THESIS_2001\Base_Line_Codes\RF2Ray  F2ray  -ASCII 
save  H:\THESIS_2001\Base_Line_Codes\3DHanning  H3D  -ASCII 

o 

*  Plot  results 

o 


figure ( 1 ) 
mesh (H3D) 
view (-37 . 5,  70) 

figure  (2 ) 

y=0 :Dz: (Nz)*Dz/2; 

plot (F, y, F2ray, y, 'o-r') 

axis (  [-30  6  0  Nz*Dz/2]) 

legend ( ' 3D  Parabolic ' , ' 2Ray ' ) 

toe 
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APPENDIX  B:  A  3D  PE  MODEL  AND  4RAY  MODEL 


Below  is  the  MatLab  program  listing  of  the  3D  parabolic  equation  model  and  the 
4Ray  model.  This  program  will  compute  the  propagation  factors  for  the  3D  PE  model 
and  the  4Ray  model  over  perfectly  flat  earth  and  constant  admittance  plane.  A  single 
absorbing  knife-edge  is  placed  between  he  transmitter  and  receiver.  The  operating 
frequency  is  1  GHz. 


o 

%  Declaring  Contants 

o 


hold  off; 
tic 
c=3e8 ; 
f=le9; 

Lamda=c/ f ; 

ko= (2*pi) /Lamda; 

S=400 ; 

Er=l  0 ; 


%  Speed  of  light 
%  Operating  frequency 
%  Wavelength 
%  Wavenumber 

%  Ground  conductivity  in  mS/m 
%  Relative  dielectric  constant 


Dz=Lamda; 

Dy=Lamda; 

sigmaz=Lamda; 


%  Delta  Z  in  meter 
%  Delta  Y  in  meter 
%  Standard  deviation  of  Z 


o 

%  Input  Parameters 

o 


Ht=60 ; 
Hk=50 ; 
Wk=4  9.8; 


%  Transmitting  antenna  height 
%  Knife  edge  height 
%  Knife  edge  width 


dl=125; 
d2=375 ; 
D=500 ; 
Dx=125 ; 


%  Distance  from  transmitter  to  knife 
%  Distance  from  knife  edge  to  receiver 
%  Range  from  the  transmitter 
%  Incremental  range  (Dx)  in  meter 


Ny=1024 ; 
Nz=l 02  4 ; 


%  Sample  size  in  the  y-direction 
%  Sample  size  in  the  z-direction 


o 

%  Define  range  in  y-direction  and  z-direction 

o 


ymax=Ny*Dy/2 ;  % 
yl=0 : Dy : ymax;  % 
y=[yl  zeros ( 1 , Ny/2-1 )] ;  % 


maximum  range  in  y-direction 
First  half  of  y  range 
y  array 


zmax  =  Nz*Dz/2; 

zl=0 : Dz : zmax; 

z=[zl  zeros ( 1 , Nz/2-1 )] ; 


%  Find  maximum  z  range  (Zmax) 

%  First  half  of  the  range  of  Zmax+1 
%  Construct  a  full  z  array 
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S^iririririritiritiri?iririririririri?*iriri?'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%  Define  Wavenumbers  in  Spatial  Domain 

S^*i(ici(irifi?ifiri(i?i(ici(*ifi?i(i?i(*if*i(irif'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 


kymax=pi/Dy; 

Dky=2  *kymax/Ny ; 
kyl=-kymax : Dky : kymax; 
ky=[kyl(:, (Ny/2 ) +1 : Ny+1 )  kyl 


%  Maximum  wavenumber  in  y-direction 
%  detla  ky 
%  Range  of  ky 
, 2 :Ny/2) ] ; 


kzmax=pi/Dz ; 

Dkz=2  *kzmax/Nz ; 
kz l=-kzmax : Dkz : kzmax; 
kz=[kzl(:, (Nz/2) +1 :Nz+l)  kzl 


%  Maximum  wavenumber  in  z-direction 
%  detla  kz 
%  Range  of  kz 
,2:Nz/2) ] ; 


o 

%  Compute  Wavenumbers  in  Frequency  Domain,  kx 

o 


Ky=meshgrid (ky, 1 :Nz) ; 
Kz=meshgrid (kz, 1 :Ny)  .  '  ; 
repeat 

kx=sqrt (koA2-Ky . A2-Kz . A2)  ; 


%  Create  a  NyxNz  matrix 
%  Create  a  NyxNz  matrix 


of  ky  row-repeat 
of  kz  column- 


%  Compute  theNy  x  Nz  kx  matrix 


clear  Ky; 
clear  Kz; 


%  Clear  matrices  Ky  and  Kz  from  memory 
%  To  make  the  algorithm  run  faster 


o 

%  Compute  reflection  coefficient 

o 


Erc=Er+i*18*S/ (f/le6)  ; 
Zs=l/sqrt (Ere) ; 

Gamma= (kz-ko*Zs) ./ (kz+ko*Zs) ; 
%gamma=-ones (Nz,Ny) ; 
gamma=meshgrid (Gamma, 1 :Ny) . ' ; 


%  Complex  dielectric 
%  Impedance 

%  Complex  reflection  coefficient 
%  Setting  gamma  =  -1 

%  Create  reflection  coefficient  matrix 
%  by  column  repeating 


S^iri?i?i?i?i?iriri?i?iririri?i?iri?'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%  Define  the  Hanning  Window 

o 


Hya= [ ] ;  % 

for  t=0:Ny/2;  % 

if  (t>=0  &  t<=3*Ny/8)  % 

hy=l ; 

elseif  (t>=3*Ny/8  &  t<=Ny/2)% 
hy= (sin (4*pi*t/Ny) ) A2; 

end 

Hya= [Hya  hy] ;  % 

end  % 

Yy=fliplr (Hya ( : , 2 :Ny/2)  )  ;  % 


An  empty  vector 
Define  Ny/2  +  1  points 
Define  first  3Ny/8  +1  points 

Define  next  Ny/8  points 

Construct  the  first  half  of  Hanning 
window  of  l+Ny/2  elements 

Flip  left 


second  element  to 
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to  right  of  the  Hanning 
%  windows  from  the 


HY= [Hya  Yy] ; 


%  Ny/2  element  for  total  of  [lx(Ny/2-l)] 
%  Hanning  window  in  y-direction 


Hzb=  [  ] ;  % 

for  t=0:Nz/2;  % 

if  ( t >— 0  &  t<=3*Nz/8)  % 

hz  =  l ; 

elseif  (t>=3*Nz/8  &  t<=Nz/2)% 
hz=(sin(4*pi*t/Nz) ) A2; 

end 

Hzb= [Hzb  hz] ;  % 

end  % 


An  empty  vector 
Define  Nz/2  +  1  points 
Define  first  3Nz/8  +1  points 

Define  next  Nz/8  points 


Construct  the  first  half  of  Hanning 
window  of  l+Nz/2  elements 


Yz=fliplr(Hzb(:,2:Nz/2) )  ; 


HZ= [Hzb  Yz] ' ; 


%  Flip  left  to  right  of  the  Hanning 
%  windows  from  the  second  element  to 
%  Nz/2  element  for  total  of  [lx(Nz/2-l)] 
%  Hanning  window  in  z-direction 


Hmy=meshgrid (HY, 1 : Nz )  ; 
Hmz=meshgrid (HZ , 1 : Ny ) 


%  Row  repeat 
%  Column  repeat 


H3D= (Hmy . *Hmz )  ; 


%  3D  Hanning  window 


o 

%  Gaussian  Current  Source 

o 


f= (1/ (sqrt (2*pi) *sigmaz) ) *exp (- (z-Ht) . A2/ (2*sigmazA2) ) ; 
f_tilda=Dz*fft (f ) ;  %  Fourier  transform  of  f(z) 

o 

%  Initial  H-field  at  x=0 

o 


g_tilda=exp (-kz . A2*sigmazA2/2) ;  %  initial  g  tilda 

hyeO_tilda=g_tilda . *cos (kz*Ht) ;  %  Initial  hye_tilda (0+, ky, kz)  field 

hyoO_tilda=-i*g_tilda . *sin (kz*Ht) ;  %  Initial  hyo_tilda (0+, ky, kz)  field 

%  Initial  even  H-field,  column  repeat 
HyeO_tilda=meshgrid (hyeO_tilda, 1 : Ny ) . 

%  Initial  odd  H-field,  column  repeat 
HyoO_tilda=meshgrid (hyoO_tilda, 1 :Ny) . 

%  Include  the  reflection  coefficient 
Hy0_tilda=0 . 5* (1+gamma) . *Hye0_tilda+0 . 5* (1-gamma) . *HyoO_tilda; 

Hy_tilda=HyO_tilda . *H3D; %  Hy_tilda (0+, ky, kz)  at  x=0 

Exp2=exp (i*kx*Dx) ;  %  Marching  range 
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o 

%  3D  Parabolic  Equation  Basic  Algorithm 

o 

C=0;  %  Setting  the  counter 


for  x=0:Dx:D-Dx; 

Hy_tilda_Dx=Hy_tilda . *Exp2 ; 

%  Hy_tilda (Dx,  ky ,  kz ) 

Hy=Dkz*Dky* (Ny*Nz*ifft2 (Hy_tilda_Dx) ) / (2*pi) A2; %  ifft2  wrt  to  ky  and 


kz 

Hy_H=Hy . *H3D;  %  Apply 

the  Hamming  window  in  spatial  domain 

C=C+Dx ; 

%  Counter 

if  (C  ==  dl)  %  When  X=  125  M  zero  out  the  field 

Hy_H ( 1 : round (Hk/Dz ) , 1 : round (Wk/Dy/ 2 ) ) =0 ; 

Hy_H ( 1 : round (Hk/Dz ) , (Ny- round (Wk/Dy/ 2 ) ) : Ny ) =0 ; 

end 


Hy  z 1 =Hy_H ( 1 : N  z / 2  + 1 ,  :) ; 

%  H-field  for  z  >  0 

HyzO=flipud (Hy_H (2:Nz/2, :) ) ; 

%  H-field  for  z  <  0 

Hye=[Hyzl;  HyzO]; 

%  Even  H-field 

Hyo=[Hyzl;  -HyzO]; 

%  Odd  H-field 

Hye_tilda= (Dz*Dy) *fft2 (Hye) ; 
Hye_tilda_H=Hye_tilda . *H3D ; 

%  Take  the  FFT  of  the  even  H~ 

%  Apply  the  Hanning  window  in  freq 

Hyo_tilda= (Dz*Dy) *fft2 (Hyo)  ; 
Hyo_tilda_H=Hyo_tilda . *H3D ; 

o, 

o 

%  Take  the  FFT  of  the  odd  H~ 

%  Apply  the  Hanning  window  in  freq 
Multiply  by  the  reflection  coefficient 

Hye_tilda_g=0 . 5* (1+gamma) . *Hye_tilda_H; 

%  Multiply  by  the  reflection  coefficient 
Hyo_tilda_g=0 . 5* (1-gamma) . *Hyo_tilda_H; 

Hy_tildal=Hye_tilda_g  +  Hyo_tilda_g; %  Hy_tilda (x, ky , kz ) 


Hy_tilda=Hy_tildal . *H3D; 

%  Apply  the  Hanning  window 

end 

o 

%  Compute  the  propagation  factor  at  the  final  range 

o 

J=2;  %  First  column  final  H-field 

rho=sqrt (DA2+ ( J)  A2  )  ;  %  Distance  from  Tx  base  of  Rx 

R=sqrt ( ( z ( : , 1 : Nz /2+1 ) -Ht ) . A2  +  rhoA2);  %  Distance  from  Tx  to  Rx 


theta=asin ( rho . /R)  ; 
phi=asin ( J/rho) ; 

g_tilda_f s=exp (- ( (cos (theta) *ko) 

%  Angle  Theta 
%  Angle  Phi 

i . A2*sigmazA2) / 2) ;  %  Free-space  g_tilda 

Hyf s=-i*ko*sin (theta) *cos (phi) . *g_tilda_fs . *exp (i*ko*R) ./ (4*pi*R) ;  % 
Free-space  H-field 
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A=Hy_H (1 : N  z / 2  + 1 , 2) ; 
F=20*logl0 (abs (A. /Hyfs ' ) ) ; 


1st  column  of  the  final  H-field 
Propagation  factor 


o, 

o 

o, 

o 


save  H:\THESIS_2001\MATLAB_CODES\PERF4Ray  F  -ASCII 

o 

%  Four-Ray  Model 

o 


%  The  radius  of  any  member  n  of  the 
perpendicular  to  the  LOS  path 

Hl=sqrt ( (Lamda*dl*d2) / (dl+d2) )  ; 

hl=Hk- (Ht*d2+z ( : , 1 :Nz/2+l) *dl) /D; 
Gamma 1=1 ; 
xl=sqrt (2) *hl/Hl; 

Cl=mfun ( ' FresnelC ' , xl ) ; 

Sl=mfun (' FresnelS ', xl )  ; 

Fnl=Cl— j*Sl; 

F1=0 . 5* (l-(l+j) *Fnl) * Gamma 1 ; 
r l=sqrt (DA2+(Ht-z(:,l:Nz/2+l) )  .  A2 )  ; 
E01=exp (- j*ko*rl) ,/rl; 

Edl=E0 1 . *F1; 

h2=Hk+ (Ht*d2-z ( : , 1 :Nz/2+l) *dl) /D; 
psiA=atan ( (Hk+Ht) /dl) ; 

%GammaA= (sin (psiA) -Zs) / (sin (psiA) +Z 
GammaA=-l ; 
x2=sqrt (2) *h2/Hl; 

C2=mfun ( ' FresnelC ' , x2 ) ; 

S2=mfun (' FresnelS x2 ) ; 

Fn2=C2- j  *S2 ; 

F2=0 . 5* (l-(l+j) *Fn2 ) *GammaA; 
r2=sqrt (DA2+(Ht+z(:,l:Nz/2+l) ) .  A2 ) ; 
E02=exp (- j *ko*r2 ) ./r2; 

Ed2=E02 . *F2 ; 

h3=Hk- (Ht*d2-z ( : , 1 :Nz/2+l) *dl) /D; 
psiB=atan ( (Hk+z (:,l:Nz/2+l))/d2); 

%  GammaB= ( sin (psiB) -Zs ) ./ (sin(psiB) 
GammaB=-l ; 
x3=sqrt (2) *h3/Hl; 

C3=mfun ( 'FresnelC ' , x3 ) ; 

S3=mfun (' FresnelS x3 )  ; 

Fn3=C3- j*S3; 

F3=0 . 5* (l-(l+j) *Fn3 ) ,*GammaB; 
r3=sqrt (DA2+(Ht  +  z(:,l:Nz/2  +  l))  .A2); 
E03=exp (- j*ko*r3) ./r3; 

Ed3=E03 . *F3 ; 

h4=Hk+ (Ht*d2+z ( : , 1 :Nz/2+l) *dl) /D; 

GammaC=GammaA*GammaB; 

x4=sqrt (2) *h4/Hl; 

C4=mf un ( 'FresnelC ' , x4 ) ; 

S4=mfun (' FresnelS ', x4 )  ; 

Fn4=C4- j  *S4 ; 


family  of  ellipsoids  in  a  plane 


%  1st  clearance  height 
%  Reflection  coefficient  for  case  1 

%  Real  part  of  the  Fresnel  Integral 
%  Imag  part  of  the  Fresnel  Integral 
%  Solution  for  the  Fresnel  Integral 
%  First  knife-edge  normalized  field 


%  2nd  clearance  height 
%  Grazing  angle  of  case  2 
s);%  Refl.  coefficient  for  case  2 


%  Real  part  of  the  Fresnel  Integral 
%  Imag  part  of  the  Fresnel  Integral 
%  Solution  for  the  Fresnel  Integral 
%  2nd  knife-edge  normalized  field 


%  3rd  clearance  height 
%  Grazing  angle  of  case  3 
+Zs);  %  Refl.  coefficient  for  case  3 

%  Real  part  of  the  Fresnel  Integral 
%  Imag  part  of  the  Fresnel  Integral 
%  Solution  for  the  Fresnel  Integral 
%  3rd  knife-edge  normalized  field 


%  4th  clearance  height 
%  Reflection  coefficient  for  case  4 

%  Real  part  of  the  Fresnel  Integral 
%  Imag  part  of  the  Fresnel  Integral 
%  Solution  for  the  Fresnel  Integral 
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F4=0 . 5* ( 1— ( 1+ j ) *Fn4 ) . *GammaC; 
r4=sqrt (DA2+(Ht-z(:,l:Nz/2+l) ) .  A2 ) ; 
E04=exp (- j *ko*r4 ) ,/r4; 

Ed4=E04 . *F4 ; 


4th  knife-edge  normalized  field 


Fknife= (Edl+Ed2+Ed3+Ed4 ) ./E01; 


Total  normalized  fields 


FknifedB=20*logl0 (abs (Fknife) ) ; 


Total  normalized  fields  in  dB 


save  H:\THESIS_2001\MATLAB_CODES\RF4Ray  FknifedB  -ASCII 


Plot  results 


figure ( 1 ) 

y=0 :Dz: (Nz)*Dz/2; 

plot (F, y, FknifedB, y, ' o-r ' ) 

axis  (  [-30  6  0  Nz*Dz/2]) 

%xlabel (' Relative  Propagation  Factor  F (dB) '), ylabel ( 'Altitude (m) ' ) 
%title(['Ny  =  'num2str(Ny)  '  Nz  =  '  num2str(Nz)  '  Dx  =  'num2str(Dx)  ' 
meters'  '  Ht  =  '  num2str(Ht)  '  meters']) 
legend ( ' Parabolic ' , ' KnifedB ' ) 
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APPENDIX  C:  A  3D  PE  MODEL  FOR  THE  MULTIPLE  SCREENS 
OF  UNIFORM  HEIGHTS  AND  EQUAL  SPACING 


Below  is  the  MatLab  program  listing  of  the  3D  parabolic  equation  model  use  to 
compute  the  propagation  factors  with  multiple  absorbing  screens  of  uniform  heights  and 


equal  spacing  over  perfectly  flat  earth  and  constant  admittance  plane.  The  operating 


frequency  is  1  GHz  and  900  MHz. 


o 

%  Declaring  Contants 

o 


hold  off; 

tic 

c=3e8 ; 

f=0 . 9e9; 

Lamda=c/f ; 

ko= (2*pi) /Lamda; 

S=400 ; 

Er=l  0 ; 


%  Speed  of  light 
%  Operating  frequency 
%  Wavelength 
%  Wavenumber 

%  Ground  conductivity  in  mS/m 
%  Relative  dielectric  constant 


Dz=Lamda; 

Dy=Lamda; 

sigmaz=Lamda; 


Delta  Z 
Delta  Y 
Current 


in  meter 
in  meter 

source  standard  deviation 


S^iritiritiritiritiritirit'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%  Input  Parameters 

o 


Ht=125 ; 


%  Transmitting  antenna  height 


Hk=20 ; 
Wk=50 ; 


%  Knife  edge  height 
%  Knife  edge  width 


NS=120 ; 

Dns=50 ; 

M=50 : Dns : Dns*NS; 


%  Number  of  screens 
%  Screen  separation 
%  Location  of  the  screens 


D=Dns* (NS+1 ) ; 
Dx=2  5 ; 


%  Range  from  the  transmitter 
%  Incremental  range  (Dx)  in  meter 


Ny=1024;  %  Sample  size  in  the  y-direction 

Nz=1024;  %  Sample  size  in  the  z-direction 

o 

%  Define  range  in  y-direction  and  z-direction 

o 


ymax=Ny*Dy/2 ;  % 
yl=0 : Dy : ymax;  % 
y=[yl  zeros ( 1 , Ny/2-1 )] ;  % 


maximum  range  in  y-direction 
First  half  of  y  range 
y  array 
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zmax  =  Nz*Dz/2;  %  Find  maximum  z  range  (Zmax) 

zl=0 : Dz : zmax;  %  First  half  of  the  range  of  Zmax+1 

z=[zl  zeros ( 1 , Nz /2-1 )] ;  %  Construct  a  full  z  array 


S^it'k-k'k'k-k-k-k-k-k-kic-k-k-k-k-k'k-kic'k-k'k-k'k'k-kic-k-k'k-k'k'k-k-k'k-k'kic-k'k-k-k-k-k-k-k-k'k-k-k'k-k-kic-k'k-kic'k-k-k-k'k'k-k'k-k'k 

o 

%  Define  Wavenumbers  in  Spatial  Domain 

S^ici(i?i?i?it*ificifi?iti?it'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

kymax=pi/Dy;  %  Maximum  wavenumber  in  y-direction 

Dky=2 *kymax/Ny ;  %  detla  ky 

kyl=-kymax : Dky : kymax;  %  Range  of  ky 

ky= [ kyl ( :  ,  (Ny/2) +1 :Ny+l)  kyl ( : , 2 : Ny/2 ) ] ; 

kzmax=pi/Dz;  %  Maximum  wavenumber  in  z-direction 

Dkz=2 *kzmax/Nz ;  %  detla  kz 

kz l=-kzmax : Dkz : kzmax;  %  Range  of  kz 

kz= [ kz 1 ( :  ,  (Nz/2) +1 :Nz  +  l)  kzl ( : , 2 : Nz/2 ) ] ; 

o 

%  Compute  Wavenumbers  in  Frequency  Domain,  kx 


Ky=meshgrid (ky,  1 

: Nz )  ; 

o, 

o 

Create  a  NyxNz  matrix 

of 

ky  row-repeat 

Kz=meshgrid (kz ,  1 
repeat 

: Ny ) . ' ; 

q. 

o 

Create  a  NyxNz  matrix 

of 

kz  column- 

kx=sqrt (koA2-Ky. 

A2-Kz . A2) ; 

o. 

o 

Compute  theNy  x  Nz  kx 

mat 

rix 

clear  Ky; 

o. 

o 

Clear  matrices  Ky  and 

Kz 

from  memory 

clear  Kz; 

o_ 

o 

To  make  the  algorithm 

run 

faster 

S^*iti?i?iriti?i?i?i?*i?i?i?*'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%  Compute  reflection  coefficient 

o 

Erc=Er+i*18*S/ (f /le6)  ;  %  Complex  dielectric 

Zs=l/sqrt (Ere) ;  %  Impedance 

Gamma= (kz-ko*Zs)  . /  (kz+ko*Zs) ;  %  Complex  reflection  coefficient 
%gamma=-ones (Nz , Ny ) ;  %  Setting  gamma  =  -1 

gamma=meshgrid (Gamma, 1 :Ny) %  Create  reflection  coefficient  matrix 

%  by  column  repeating 

o 

%  Define  the  Hanning  Window 


An  empty  vector 
Define  Ny/2  +  1  points 
Define  first  3Ny/8  +1  points 


Hya=  [  ] ; 
for  t=0:Ny/2; 

if  ( t >— 0  &  t<=3*Ny/8) 
hy=l  ; 

elseif  (t>=3*Ny/8  &  t<=Ny/2)%  Define  next  Ny/8  points 
hy= (sin (4*pi*t/Ny) ) A2; 

end 

Hya= [Hya  hy] ; 

end 


%  Construct  the  first  half  of  Hanning 
%  window  of  l+Ny/2  elements 
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Yy=f liplr (Hya ( : , 2 : Ny/2  )  )  ; 
second  element  to 
HY= [Hya  Yy] ; 


%  Flip  left  to  right  of  the  Hanning 

%  windows  from  the 

%  Ny/2  element  for  total  of  [lx(Ny/2-l)] 
%  Hanning  window  in  y-direction 


Hzb=  [  ] ;  % 

for  t=0:Nz/2;  % 

if  ( t >— 0  &  t<=3*Nz/8)  % 

hz  =  l ; 

elseif  (t>=3*Nz/8  &  t<=Nz/2)% 
hz=(sin(4*pi*t/Nz) ) A2; 

end 

Hzb= [Hzb  hz ] ;  % 

end  % 


An  empty  vector 
Define  Nz/2  +  1  points 
Define  first  3Nz/8  +1  points 

Define  next  Nz/8  points 

Construct  the  first  half  of  Hanning 
window  of  l+Nz/2  elements 


Yz=fliplr(Hzb(:,2:Nz/2) )  ; 


HZ= [Hzb  Yz] ' ; 


%  Flip  left  to  right  of  the  Hanning 
%  windows  from  the  second  element  to 
%  Nz/2  element  for  total  of  [lx(Nz/2-l)] 
%  Hanning  window  in  z-direction 


Hmy=meshgrid (HY, 1 : Nz )  ; 
Hmz=meshgrid (HZ , 1 : Ny )  '  ; 


%  Row  repeat 
%  Column  repeat 


H3D= (Hmy . *Hmz ) ; 


%  3D  Hanning  window 


o 

%  Gaussian  Current  Source 

S^ici(*i?i?iti?ificifi?i?iciti?it'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 


f= (1/ (sqrt (2*pi) *sigmaz) ) *exp (- (z-Ht) . A2/ (2*sigmazA2) ) ; 
f_tilda=Dz*f ft (f ) ;  %  Fourier  transform  of  f(z) 

o 

%  Initial  H-field  at  x=0 

o 


g_tilda=exp (-kz . A2*sigmazA2/2) ;  %  initial  g  tilda 

hyeO_tilda=g_tilda . *cos (kz*Ht) ;  %  Initial  hye_tilda (0+, ky, kz)  field 

hyoO_tilda=-i*g_tilda . *sin (kz*Ht) ;  %  Initial  hyo_tilda (0+, ky, kz)  field 

%  Initial  even  H-field,  column  repeat 
HyeO_tilda=meshgrid (hyeO_tilda, 1 :Ny) . 

%  Initial  odd  H-field,  column  repeat 
HyoO_tilda=meshgrid (hyoO_tilda, 1 :Ny) . 

%  Include  the  reflection  coefficient 
Hy0_tilda=0 . 5* (1+gamma) . *Hye0_tilda+0 . 5* (1-gamma) . *HyoO_tilda; 

Hy_tilda=HyO_tilda . *H3D; %  Hy_tilda (0+, ky, kz)  at  x=0 

Exp2=exp (i*kx*Dx) ;  %  Marching  range 
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o 


%  3D  Parabolic  Equation  Basic  Algorithm 

o 

C=0;  %  Setting  the  counter 


for  x=0:Dx:D-Dx; 

Hy_tilda_Dx=Hy_tilda . *Exp2 ; 

%  Hy_tilda (Dx,  ky ,  kz ) 

%  ifft2  wrt  to  ky  and  kz 

Hy=Dkz*Dky* (Ny*Nz*ifft2 (Hy_tilda_Dx) ) / (2*pi) A2; 

Hy_H=Hy . *H3D;  %  Apply  the  Hamming  window  in  spatial  domain 

C=C+Dx;  %  Counter 

while  M(M==C) 

Hy_H  ( 1 : round (Hk/Dz ) , 1 : round (Wk/Dy/ 2 ) ) =0  ; 

Hy_H ( 1 : round (Hk/Dz ) , (Ny-round (Wk/Dy/ 2) ) : Ny ) =0 ; 
break 

end 


Hyzl=Hy_H (1 :Nz/2+l, :) ; 

%  H-field  for  z  >  0 

HyzO=flipud (Hy_H (2:Nz/2,  :) )  ; 

%  H-field  for  z  <  0 

Hye=[Hyzl;  HyzO]; 

%  Even  H-field 

Hyo=[Hyzl;  -HyzO]; 

%  Odd  H-field 

Hye_tilda= (Dz*Dy) *fft2 (Hye)  ; 
Hye_tilda_H=Hye_tilda . *H3D ; 

%  Take  the  FFT  of  the  even  H~ 

%  Apply  the  Hanning  window  in  freq 

Hyo_tilda= (Dz*Dy) *fft2 (Hyo) ; 
Hyo_tilda_H=Hyo_tilda . *H3D ; 

o. 

o 

%  Take  the  FFT  of  the  odd  H~ 

%  Apply  the  Hanning  window  in  freq 
Multiply  by  the  reflection  coefficient 

Hye_tilda_g=0 . 5* (1+gamma) . *Hye_tilda_H; 

%  Multiply  by  the  reflection  coefficient 
Hyo_tilda_g=0 . 5* ( 1-gamma) . *Hyo_tilda_H; 

Hy_tildal=Hye_tilda_g  +  Hyo_tilda_g; %  Hy_tilda (x, ky , kz ) 
Hy_tilda=Hy_tildal . *H3D;  %  Apply  the  Hanning  window 

end 

o 

%  Compute  The  Propagation  Factor 

^A’A’A’A’A,A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A,A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’A’,A'A'A'A'A'A 

o 


II 

%  First  column  final  H-field 

rho=sqrt (DA2+ (J) A2) ; 

%  Distance  from  Tx  to  base  Rx 

R=sqrt ( (z ( : , 1 :Nz/2+l) -Ht) . A2  +  rhoA2);  %  Distance  from  Tx  to  Rx 


theta=asin ( rho . /R) ; 
phi=asin ( J/ rho) ; 

g_tilda_f s=exp (- ( (cos (theta) *ko) 

%  Angle  Theta 
%  Angle  Phi 

i . A2*sigmazA2) /2) ;  %  Free-space  g_tilda 

%  Free-space  H-field 

Hyf s=-i*ko*sin (theta) *cos (phi)  . *g_tilda_fs . *exp (i*ko*R)  ./ (4*pi*R)  ; 
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A=Hy_H (1 :Nz/2+l, 2) ;  %  z  >  0  of  the  first  column  of  the  final  H-field 
F= (abs (A. /Hyfs . ' ) ) ;  %  Propagation  factor  relative  to  the  free-space 

save  H:\THESIS_2001\MATLAB_CODES\propfact  F  -ASCII 

o 

%  Plot  results 

S^icititititititifititititic'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

figure ( 1 ) 

%Z=linspace (0, 3*Nz/8*Dz,  513)  ; 

Z= (0 :Dz :Nz/2*Dz)  '  ; 

Y= (Z-Hk) /Lamda; 

save  H:\THESIS_2001\MATLAB_CODES\recantenna  Y  -ASCII 
plot (Y, F) 

axis (  [-30  60  0  1.8  ]) 

ylabel (' Relative  Propagation  Factor,  F '), xlabel (' Receiving  Antenna  in 
Wavelegth ' ) 

title (['Ny  =  'num2str(Ny)  '  Nz  =  '  num2str(Nz)  '  Dx  =  'num2str(Dx)  ' 

meters'  '  Ht  =  '  num2str(Ht)  '  meters']) 

grid 

toe  %  Stop  stopwatch 
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APPENDIX  D:  A  3D  PE  MODEL  FOR  A  SINGLE  ROUND  HILL 


Below  is  the  MatLab  program  listing  of  the  3D  parabolic  equation  model  use  to 
compute  the  field  strengths  for  a  single  round  hill  with  multiple  absorbing  screens  of 
uniform  heights  and  equal  spacing.  These  screens  represent  a  row  of  buildings  or  houses 
in  residential  areas  of  a  well  built-up  city.  The  operating  frequency  is  900  MHz. 

o 

%  Declaring  Contants 

o 


hold  off; 
tic 
c=3e8 ; 
f=0 . 9e9; 

Lamda=c/ f ; 

ko= (2*pi) /Lamda; 

S=400 ; 

Er=l  0 ; 


%  Speed  of  light 
%  Operating  frequency 
%  Wavelength 
%  Wavenumber 

%  Ground  conductivity  in  mS/m 
%  Relative  dielectric  constant 


Dz=Lamda; 

Dy=Lamda; 

sigmaz=Lamda; 


%  Delta  Z  in  meter 
%  Delta  Y  in  meter 

%  Current  source  standard  deviation 


%  Input  Parameters 


Ht=57 ; 
Hbs=7 ; 


%  Transmitting  antenna  height 
%  Receiving  antenna  height  above  KE ' s 


Wk=100 ; 
Dns=50 ; 


%  Knife  edge  width  of  25  meter 
%  Screen  separation 


Dist=10000 ;  % 
xs=-1000 : Dns : 1000 ;  % 
ys=sqrt (le4A2-xs.^2)-(9.95e3) ;% 
xsl=1000+Dns : Dns : Dist+Dns;  % 
ysl  =  zeros (1,  (length (xsl) ) ) ;  % 
xs2= [xs  xsl ] ; 

Hk= [ys  ysl ] ;  % 
Hk (1) =0;  % 
Hk ( length (xs ) ) =0 ; 


Final  range 
Horizontal  hill  width 
Equation  for  circle 

Equal  heigth  and  equal  distance  screen 
Equal  heigth  screen  height 

Knife  edge  height  of  meter 

When  screen  height  less  than  zeros 


NS=length (xs2 ) ; 
M=Dns :Dns:Dns* (NS) ; 


%  Number  of  screens 
%  Location  of  the  screens 


D=Dns* (NS) ; 
Dx=50 ; 


%  Range  from  the  transmitter 
%  Incremental  range  (Dx)  in  meter 


Ny=1024 ; 
Nz  =  l 02  4 ; 


%  Sample  size  in  the  y-direction 
%  Sample  size  in  the  z-direction 
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9?iri?iriti?i?iriri?i?iri?iri?iritiri?iritiri?iriririr'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 


%  Define  range  in  y-direction  and  z-direction 

S^*itiriti?iti?iti?i?i?i?i?iti?iti?ifiri?*i(iriti?i('k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 


ymax=Ny*Dy/2 ;  % 
yl=0 : Dy : ymax;  % 
y=[yl  zeros ( 1 , Ny/2-1 )] ;  % 


maximum  range  in  y-direction 
First  half  of  y  range 
y  array 


zmax  =  Nz*Dz/2; 

zl=0 : Dz : zmax; 

z— [ z 1  zeros (l,Nz/2-l) ] ; 


%  Find  maximum  z  range  (Zmax) 

%  First  half  of  the  range  of  Zmax+1 
%  Construct  a  full  z  array 


o 

%  Define  Wavenumbers  in  Spatial  Domain 

o 


kymax=pi/Dy; 

Dky=2  *kymax/Ny ; 
kyl=-kymax : Dky : kymax; 
ky=[kyl(:, (Ny/2 ) +1 : Ny+1 )  kyl ( 

kzmax=pi/Dz ; 

Dkz=2  *kzmax/Nz ; 

kz l=-kzmax : Dkz : kzmax; 

kz= [kzl ( : , (Nz/2) +1 :Nz+l)  kzl ( 


%  Maximum  wavenumber  in  y-direction 
%  detla  ky 
%  Range  of  ky 
, 2 :Ny/2) ] ; 

%  Maximum  wavenumber  in  z-direction 
%  detla  kz 
%  Range  of  kz 
, 2 :Nz/2) ] ; 


o 

%  Compute  Wavenumbers  in  Frequency  Domain,  kx 

o 


Ky=meshgrid (ky, 1 :Nz)  ; 
Kz=meshgrid (kz, 1 :Ny)  .  '  ; 
repeat 

kx=sqrt (koA2-Ky . A2-Kz . A2)  ; 


%  Create 
%  Create 


a  NyxNz  matrix 
a  NyxNz  matrix 


of  ky  row-repeat 
of  kz  column- 


%  Compute  theNy  x  Nz  kx  matrix 


clear  Ky; 
clear  Kz; 


%  Clear  matrices  Ky  and  Kz  from  memory 
%  To  make  the  algorithm  run  faster 


o 

%  Compute  reflection  coefficient 

o 


Erc=Er+i*18*S/ (f/le6)  ; 
Zs=l/sqrt (Ere)  ; 

Gamma= (kz-ko*Zs)  ./  (kz+ko*Zs) ; 
%gamma=-ones (Nz,Ny)  ; 
gamma=meshgrid (Gamma, 1 :Ny) . ' ; 


%  Complex  dielectric 
%  Impedance 

%  Complex  reflection  coefficient 
%  Setting  gamma  =  -1 

%  Create  reflection  coefficient  matrix 
%  by  column  repeating 


o 

%  Define  the  Hanning  Window 

o 


Hya=  [  ]  ; 

for  t=0:Ny/2; 


%  An  empty  vector 
%  Define  Ny/2  +  1  points 
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if  ( t >— 0  &  t<=3*Ny/8) 
hy=l ; 

elseif  (t>=3*Ny/8  &  t<=Ny/2 
hy= (sin (4*pi*t/Ny) )  A2; 

end 

Hya= [Hya  hy] ; 

end 


%  Define  first  3Ny/8  +1  points 
%  Define  next  Ny/8  points 

%  Construct  the  first  half  of  Hanning 
%  window  of  l+Ny/2  elements 


Yy=fliplr (Hya (  :  ,  2 :Ny/2)  )  ; 
second  element  to 
HY= [Hya  Yy] ; 


%  Flip  left  to  right  of  the  Hanning 
%  windows  from  the 

%  Ny/2  element  for  total  of  [lx(Ny/2-l)] 
%  Hanning  window  in  y-direction 


Hzb=  [  ] ;  % 

for  t=0:Nz/2;  % 

if  (t>=0  &  t<=3*Nz/8)  % 

hz  =  l ; 

elseif  (t>=3*Nz/8  &  t<=Nz/2)% 
hz=(sin(4*pi*t/Nz)  )  A2; 

end 

Hzb= [Hzb  hz ] ;  % 

end  % 


An  empty  vector 
Define  Nz/2  +  1  points 
Define  first  3Nz/8  +1  points 

Define  next  Nz/8  points 

Construct  the  first  half  of  Hanning 
window  of  l+Nz/2  elements 


Yz=fliplr(Hzb(:,2:Nz/2)  )  ; 


HZ= [Hzb  Yz] ' ; 


%  Flip  left  to  right  of  the  Hanning 
%  windows  from  the  second  element  to 
%  Nz/2  element  for  total  of  [lx(Nz/2-l)] 
%  Hanning  window  in  z-direction 


Hmy=meshgrid (HY,  1 : Nz )  ; 
Hmz=meshgrid (HZ , 1 : Ny ) 


%  Row  repeat 
%  Column  repeat 


H3D= (Hmy . *Hmz ) ; 


%  3D  Hanning  window 


o 

%  Gaussian  Current  Source 

O 


f= (1/ (sqrt (2*pi) *sigmaz) ) *exp (- (z-Ht) . A2/ (2*sigmazA2) ) ; 
f_tilda=Dz*fft (f ) ;  %  Fourier  transform  of  f(z) 

o 

%  Initial  H-field  at  x=0 

o 


g_tilda=exp (-kz . A2*sigmazA2/2) ;  %  initial  g  tilda 

hyeO_tilda=g_tilda . *cos (kz*Ht) ;  %  Initial  hye_tilda (0+, ky, kz)  field 

hyoO_tilda=-i*g_tilda . *sin (kz*Ht) ;  %  Initial  hyo_tilda (0+, ky, kz)  field 

%  Initial  even  H-field,  column  repeat 
HyeO_tilda=meshgrid (hyeO_tilda, 1 :Ny) . 

%  Initial  odd  H-field,  column  repeat 
HyoO_tilda=meshgrid (hyoO_tilda, 1 :Ny) . 

%  Include  the  reflection  coefficient 
Hy0_tilda=0 . 5* (1+gamma) . *Hye0_tilda+0 . 5* (1-gamma) . *HyoO_tilda; 

Hy_tilda=HyO_tilda . *H3D; %  Hy_tilda (0+, ky, kz)  at  x=0 
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Exp2=exp (i*kx*Dx)  ; 


Marching  range 


o, 

o 


o 

%  3D  Parabolic  Equation  Algorithm 

o 

C=0;  %  Setting  the  range  counter 

c=0;  %  Setting  the  index  counter 


for  x=0:Dx:D-Dx; 
c=c+l ; 


%  Define  the  total  steps  based  on  Dx  and  D 
%  Starting  the  index  counter 


Hy_tilda_Dx=Hy_tilda . *Exp2 ;  %  Hy_tilda (Dx, ky ,  kz ) 

%  ifft2  wrt  to  ky  and  kz 
Hy=Dkz*Dky* (Ny*Nz*ifft2 (Hy_tilda_Dx) ) / (2*pi) A2; 


Hy_H=Hy . *H3D;  %  Apply  the  Hamming  window  in  spatial  domain 

o 

%  Insert  Screens  (Buildings) 

o 

C=C+Dx;  %  Range  Counter 

while  M(M==C)  %  If  location  of  screen  =  marching  distance 

Hy_H  ( 1 : round ( (Hk (c) +Hbs ) /Dz ) , 1 : round (Wk/Dy/ 2 ) ) =0  ; 

Hy_H ( 1 : round ( (Hk (c) +Hbs ) /Dz ) , (Ny-round (Wk/Dy/ 2 ) ) : Ny ) =0 ; 
break 

end 

o 

%  Field  Strength  Convert  from  3D  to  2D 

S^i?i?i?*it*i?iciti?it*i?i?it*i?irif*'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%  Measure  the  field  strength  at  the  top  of  each  screen 

o 

J=4;  %  First  column  final  H-field 

rho=sqrt (CA2+ ( J) A2 )  ;  %  Distance  from  transmitter  to  the  receiver 

R=sqrt ( ( (Hk (c) +Hbs) -Ht) A2  +  rhoA2); 

%  Distance  from  transmitter  to  receiver 

theta=asin (rho/R)  ;  %  Angle  Theta  measured  relative  to  the  z-axis 

phi=asin ( J/rho) ;  %  Angle  Phi  measured  relative  to  the  x-axis 

%  Free-space  g_tilda 

g_tilda_f s=exp (- ( (cos (theta) *ko) . A2*sigmazA2) /2) ; 

%  Free-space  H-field 

Hyfs=-i*ko*sin (theta) *cos (phi) . *g_tilda_f s . *exp (i*ko*R) ./ (4*pi*R) ; 

%  z  >  0  of  the  first  column  of  the  final  H-field 
A=Hy_H (round ( (Hk (c) +Hbs) /Dz) +1,4)  ; 

Frel (c) =20*logl0 (abs (A/Hyfs) ) ;  %  Relative  Prop  Factor 

%  Convert  propagation  factor  (F)  to  E-field  by  multiplying  it  by 
%  1/sqrt (ko*R) 

Field (c) =20*logl0 (abs (A/Hyfs) /sqrt (ko*R) ) ;  %  E-Field  strength 

FFreeSpace (c) =20*logl0 (1/sqrt (ko*R) ) ;  %  2D  free-space  pathloss 
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S^itititirirititiritir'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 


Hy z 1 =Hy_H ( 1 : N z / 2  + 1 ,  :)  ; 

HyzO=flipud (Hy_H (2 :Nz/2,  :)  )  ; 

Hye=[Hyzl;  HyzO]; 

Hyo=[Hyzl;  -HyzO]; 

Hye_tilda= (Dz*Dy) *fft2 (Hye)  ; 
Hye_tilda_H=Hye_tilda . *H3D; 

Hyo_tilda= (Dz*Dy) *fft2 (Hyo)  ; 
Hyo_tilda_H=Hyo_tilda . *H3D; 

%  Apply  reflection 
Hye_tilda_g=0 . 5* ( 1+gamma) . *Hye_ti 


H-field  for  z  >  0 
H-field  for  z  <  0 
Even  H-field 
Odd  H-field 

Take  the  FFT  of  the  even  H-tilda 
the  Hanning  window  in  freq  domain 

Take  the  FFT  of  the  odd  H-tilda 
The  Hanning  window  in  freq  domain 
oefficient  to  the  even  field 
.a_H; 


%  Apply  reflection  coefficiet  to  the  odd  field 
Hyo_tilda_g=0 . 5* ( 1-gamma) . *Hyo_tilda_H; 


Hy_tildal=Hye_tilda_g  +  Hyo_tilda_g;  %  Hy_tilda (x, ky , kz ) 


Hy_tilda=Hy_tildal . *H3D;  %  The  Hanning  window  in  freq  domain 


end 

save  H:\THESIS_2001\MATLAB_CODES\FieldStl00  Field  -ASCII 
save  H:\THESIS_2001\MATLAB_CODES\Free_Space  FFreeSpace  -ASCII 

o 

%  Plot  Results 

o 

figure  ( 1 ) 

stem (xs2 , Hk+Hbs ) , grid 
xlabel (' Screen  Placement  (m) ') 
ylabel ( ' Screen  Height  (m) ' ) 

figure (2 ) 

Z= (-1000 :Dns :Dist+Dns) '; 

plot ( Z , Frel ) ,grid 

axis (  [-1000  5000  -80  10] ) 

ylabel (' Relative  Field  Strength  (dB) ') 

xlabel (' Screen  Position  (m) ') 

title ([ '  Wk  =  '  num2str(Wk)  '  meters'  '  Dx  =  'num2str (Dx)  '  meters'  ' 
Ht  =  '  num2str(Ht)  '  meters']) 

figure ( 3 ) 

Z= (-1000 :Dns :Dist+Dns) '; 

save  H:\THESIS_2001\MATLAB_CODES\RHPlacement  Z  -ASCII 

plot (Z, Field, Z, FFreeSpace) , grid 

axis ( [-1000  5000  -120  -30]) 

ylabel (' Field  Strength  (dB)  ') 

xlabel (' Screen  Position  (m) ') 

title ([ '  Wk  =  '  num2str(Wk)  '  meters'  '  Dx  =  'num2str (Dx)  '  meters']) 

%  Stop  stopwatch 
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toe 
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APPENDIX  E:  A  3D  PE  MODEL  FOR  TWO  HILLS  OF 

SINUSOIDAL  SHAPE 


Below  is  the  MatLab  program  listing  of  the  3D  parabolic  equation  model  use  to 
compute  the  field  strengths  for  two  sinusoidal  hills  with  multiple  absorbing  screens  of 
uniform  heights  and  equal  spacing.  These  screens  represent  a  row  of  buildings  or  houses 
in  residential  areas  of  a  well  built-up  city.  The  operating  frequency  is  900  MHz. 

o 

%  Declaring  Contants 

o 


hold  off; 

tic 

c=3e8 ; 

f=0 . 9e9; 

Lamda=c/f ; 

ko= (2*pi) /Lamda; 

S=400 ; 

Er=l  0 ; 


%  Speed  of  light 
%  Operating  frequency 
%  Wavelength 
%  Wavenumber 

%  Ground  conductivity  in  mS/m 
%  Relative  dielectric  constant 


Dz=Lamda; 

Dy=Lamda; 

sigmaz=Lamda; 


Delta  Z 
Delta  Y 
Current 


in  meter 
in  meter 

source  standard  deviation 


S^iritiritiritiritiritirit'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%  Input  Parameters 

o 


Ht=57 ; 
Hbs=7 ; 


%  Transmitting  antenna  height 
%  Receiving  antenna  height  above  KE ' s 


Wk=100 ;  % 

Dns=50;  % 

xs=-1500 : Dns : 4500+Dns;  % 

Hk=Hbs+25+25*cos (2*pi*xs/3000)  ; 


Knife  edge  width  of  25  meter 
Screen  separation 

Horizontal  hill  width 

%  Equal  space  screen  height 


NS=length (xs )  ; 

M=Dns :Dns :Dns* (NS)  ; 


%  Number  of  screens 
%  Location  of  the  screens 


D=Dns* (NS) ; 
Dx=50 ; 


%  Range  from  the  transmitter 
%  Incremental  range  (Dx)  in  meter 


Ny=1024 ; 
Nz=l 02  4 ; 


Sample  size 
Sample  size 


in  the  y-direction 
in  the  z-direction 


o 

%  Define  range  in  y-direction  and  z-direction 


%  maximum  range  in  y-direction 
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ymax=Ny*Dy/ 2 ; 


yl=0 : Dy : ymax; 

y=[yl  zeros ( 1 , Ny/2-1 )] ; 


%  First  half  of  y  range 
%  y  array 


zmax  =  Nz*Dz/2; 

zl=0 : Dz : zmax; 

z— [ z 1  zeros (l,Nz/2-l) ] ; 


%  Find  maximum  z  range  (Zmax) 

%  First  half  of  the  range  of  Zmax+1 
%  Construct  a  full  z  array 


Q'^-k'k-k-k-k-k-k-k-k'k-k-k-k-k-k-k'k'k-k-k-k-k-k-k^-k-k-k-k-k-k-k^-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k'k-k-k-k'k-k'k-k-k^ 

o 

%  Define  Wavenumbers  in  Spatial  Domain 


kymax=pi/Dy; 

%  Maximum  wavenumber 

in  y-direction 

Dky=2  *kymax/Ny ; 

%  detla  ky 

kyl=-kymax : Dky : kymax; 

%  Range  of  ky 

ky=[kyl(:, (Ny/2 ) +1 : Ny+1 ) 

kyl ( : , 2 :Ny/2) ] ; 

kzmax=pi/Dz ; 

%  Maximum  wavenumber 

in  z-direction 

Dkz=2  *kzmax/Nz ; 

%  detla  kz 

kz l=-kzmax : Dkz : kzmax; 

%  Range  of  kz 

kz=[kzl(:, (Nz/2) +1 :Nz+l) 

kzl  (  : , 2 :Nz/2) ] ; 

o 

%  Compute  Wavenumbers  in  Frequency  Domain,  kx 

S^*i(i?i(iri(i?i(irifi?i(*i(i?i(*ifi?i(icifi?i('k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 


Ky=meshgrid (ky, 1 :Nz)  ; 
Kz=meshgrid (kz, 1 :Ny)  .  '  ; 
repeat 

kx=sqrt (koA2-Ky . A2-Kz . A2)  ; 

clear  Ky; 
clear  Kz; 


Create  a  NyxNz  matrix 
Create  a  NyxNz  matrix 


of  ky  row-repeat 
of  kz  column- 


Compute  theNy  x  Nz  kx  matrix 


Clear  matrices  Ky  and  Kz  from  memory 
To  make  the  algorithm  run  faster 


o 

%  Compute  reflection  coefficient 

S^icititifitititificifititititiritititit'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

Erc=Er+i*18*S/ (f /le6) ;  %  Complex  dielectric 

Zs=l/sqrt (Ere) ;  %  Impedance 

Gamma= (kz-ko*Zs)  . / (kz+ko*Zs)  ;  %  Complex  reflection  coefficient 
%gamma=-ones (Nz , Ny ) ;  %  Setting  gamma  =  -1 

gamma=meshgrid (Gamma, 1 :Ny) .' ;  %  Create  reflection  coefficient  matrix 

%  by  column  repeating 

o 

%  Define  the  Hanning  Window 

o 

Hya= [ ] ;  %  An  empty  vector 

for  t=0:Ny/2;  %  Define  Ny/2  +  1  points 

if  (t>=0  &  t<=3*Ny/8)  %  Define  first  3Ny/8  +1  points 

hy=l ; 

elseif  (t>=3*Ny/8  &  t<=Ny/2)%  Define  next  Ny/8  points 
hy= (sin(4*pi*t/Ny) ) A2; 

end 
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Hya= [Hya  hy] ; 

end 


%  Construct  the  first  half  of  Hanning 
%  window  of  l+Ny/2  elements 


Yy=fliplr (Hya ( : , 2 :Ny/2)  )  ; 
second  element  to 
HY= [Hya  Yy] ; 


%  Flip  left  to  right  of  the  Hanning 
%  windows  from  the 

%  Ny/2  element  for  total  of  [lx(Ny/2-l)] 
%  Hanning  window  in  y-direction 


Hzb=  [  ] ;  % 

for  t=0:Nz/2;  % 

if  (t>=0  &  t<=3*Nz/8)  % 

hz  =  l ; 

elseif  (t>=3*Nz/8  &  t<=Nz/2)% 
hz=(sin(4*pi*t/Nz)  )  A2; 

end 

Hzb= [Hzb  hz ] ;  % 

end  % 


An  empty  vector 
Define  Nz/2  +  1  points 
Define  first  3Nz/8  +1  points 

Define  next  Nz/8  points 

Construct  the  first  half  of  Hanning 
window  of  l+Nz/2  elements 


Yz=fliplr(Hzb(:,2:Nz/2) )  ; 


HZ= [Hzb  Yz] ' ; 


%  Flip  left  to  right  of  the  Hanning 
%  windows  from  the  second  element  to 
%  Nz/2  element  for  total  of  [lx(Nz/2-l)] 
%  Hanning  window  in  z-direction 


Hmy=meshgrid (HY,  1 : Nz )  ; 
Hmz=meshgrid (HZ , 1 : Ny ) ' ; 


%  Row  repeat 
%  Column  repeat 


H3D= (Hmy . *Hmz ) ; 


%  3D  Hanning  window 


o 

%  Gaussian  Current  Source 

o 


f= (1/ (sqrt (2*pi) *sigmaz) ) *exp (- (z-Ht) . A2/ (2*sigmazA2) ) ; 
f_tilda=Dz*fft (f ) ;  %  Fourier  transform  of  f(z) 

S^icitititicitititicifititititit'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%  Initial  H-field  at  x=0 

o 


g_tilda=exp  (-kz  .  ,'2*sigmazA2/2)  ;  %  initial  g  tilda 

hyeO_tilda=g_tilda . *cos (kz*Ht) ;  %  Initial  hye_tilda (0+, ky, kz)  field 

hyoO_tilda=-i*g_tilda . *sin (kz*Ht) ;  %  Initial  hyo_tilda (0+, ky, kz)  field 

%  Initial  even  H-field,  column  repeat 
HyeO_tilda=meshgrid (hyeO_tilda, 1 :Ny) . 

%  Initial  odd  H-field,  column  repeat 
HyoO_tilda=meshgrid (hyoO_tilda, 1 :Ny) . 

%  Include  the  reflection  coefficient 
Hy0_tilda=0 . 5* (1+gamma) . *Hye0_tilda+0 . 5* (1-gamma) . *HyoO_tilda; 

Hy_tilda=HyO_tilda . *H3D; %  Hy_tilda (0+, ky, kz)  at  x=0 

Exp2=exp (i*kx*Dx) ;  %  Marching  range 
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o 

%  3D  PE  Basic  Algorithm 

o 


C=0  ; 
c=0  ; 


%  Setting  the  range  counter 
%  Setting  the  index  counter 


for  x=0:Dx:D-Dx; 
D 

c=c+l ; 


%  Define  the  total  steps  based  on  Dx  and 
%  Starting  the  index  counter 


Hy_tilda_Dx=Hy_tilda . *Exp2 ; %  Hy_tilda (Dx, ky , kz ) 

%  ifft2  wrt  to  ky  and  kz 
Hy=Dkz*Dky* (Ny*Nz*ifft2 (Hy_tilda_Dx) ) / (2*pi) A2; 


Hy_H=Hy . *H3D; 


%  Apply  the  Hamming  window  in  spatial 


domain 


o 

%  Insert  Screens  (Buildings) 


C=C+Dx ; 


%  Range  Counter 


while  M(M==C)  %  If  location  of  screen  =  marching  distance 

%  only  work  if  Dx  =  Dns;  otherwise  rewrite 
Hy_H  ( 1 : round ((Hk(c))/Dz),l: round (Wk/Dy/ 2 ) ) =0  ; 

Hy_H ( 1 : round ( (Hk (c) ) /Dz ) , (Ny-round (Wk/Dy/ 2 ) ) : Ny ) =0 ; 
break 

end 

o 

%  Field  Strength  Convert  from  3D  to  2D 


J=4; 


%  First  column  final  H-field 
%  Distance  from  Tx  to  Base  Rx 


rho=sqrt (CA2  + (J) A2) ; 

R=sqrt ( ( (Hk (c) ) -Ht) A2  +  rhoA2);  %  Distance  from  Tx  to  Rx 
theta=asin ( rho/R) ;  %  Angle  Theta 

phi=asin ( J/rho) ;  %  Angle  Phi 

%  Free-space  g_tilda 

g_tilda_f s=exp (- ( (cos (theta) *ko) . A2*sigmazA2) / 2)  ; 

%  Free-space  H-field 

Hyfs=-i*ko*sin (theta) *cos (phi)  . *g_tilda_fs . *exp (i*ko*R)  ./ (4*pi*R)  ; 

%  z  >  0  of  the  first  column  of  the  final  H-field 
A=Hy_H (round ((Hk(c))/Dz)+1,4); 


Frel (c) =20*logl0 (abs (A/Hyfs) )  ; 
Field(c)=20*logl0 (abs (A/Hyfs) /sqrt (ko*R) )  ; 

FFreeSpace (c) =20*logl0 (abs (Hyfs) /sqrt (ko*R) )  ; 


%  Field  strength 
%  Field  strength 


SL;  "k  ★  ★  ★  -k  ★  -k  ★  -k  -k  ★  ★  ★  -k  ★  ★  ★  -k  ★  ★  ★  -k  ★  ■*■  ★  -k  ★  ★  ★  ★  ★  ★  ★  ~k  ★  ■*■  ★  ~k  ★  ★  ★  ~k  ★  ★  ★  ■*■  ★  -k  ★  ★  -k  ★  ★  ★  ■*■  ★  -k  ★  ★  ★  ★  ★  ★ 

%  H-field  for  z  >  0 


Hy z 1 =Hy_H ( 1 : N z / 2  + 1 ,  :) ; 
HyzO=flipud (Hy_H (2:Nz/2, :) ) ; 


%  H-field  for  z  <  0 

%  Even  H-field 
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Hye=[Hyzl;  HyzO] 


Hyo=[Hyzl;  -HyzO]; 

Hye_tilda= (Dz*Dy) *fft2 (Hye)  ; 
Hye_tilda_H=Hye_tilda . *H3D; 

Hyo_tilda= (Dz*Dy) *fft2 (Hyo)  ; 
Hyo_tilda_H=Hyo_tilda . *H3D; 

%  Apply  reflect! 
Hye_tilda_g=0 . 5* ( 1+gamma) . *Hye_ 

%  Apply  reflecti 
Hyo_tilda_g=0 . 5* ( 1-gamma) . *Hyo_ 


%  Odd  H-field 

%  Take  the  FFT  of  the  even  H-tilda 
%  the  Hanning  window  in  freq  domain 

%  Take  the  FFT  of  the  odd  H-tilda 
%  The  Hanning  window  in  freq  domain 
n  coefficient  to  the  even  field 
ilda_H; 

n  coefficiet  to  the  odd  field 
ilda_H; 


Hy_tildal=Hye_tilda_g  +  Hyo_tilda_g;  %  Hy_tilda (x,  ky ,  kz ) 


Hy_tilda=Hy_tildal . *H3D;  %  The  Hanning  window  in  freq  domain 


end 

save  H:\THESIS_2001\Base_Line_Codes\TwoHFieldl00  Field  -ASCII 


S^iritititiririrititit'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%  Plot  Results 

o 

figure  ( 1 ) 

stem (xs, Hk) , grid 

axis  (  [-2000  5000  0  60] ) 

xlabel (' Screen  Placement  (m) ') 

ylabel ( ' Screen  Height  (m) ' ) 

figure (2 ) 

Z= (-1500 :Dns : 4500+Dns) '; 

plot ( Z , Frel ) ,grid 

axis ( [-2000  5000  -80  10]  ) 

ylabel (' Relative  Field  Strength  (dB) ') 

xlabel (' Screen  Position  (m) ') 

title ([ '  Wk  =  '  num2str(Wk)  '  meters'  '  Dx  =  'num2str (Dx)  '  meters'  ' 

Ht  =  '  num2str(Ht)  '  meters']) 

figure ( 3 ) 

Z=  (-1500 :Dns : 4500+Dns)  '; 
plot (Z, Field) ,grid 
axis ( [-2000  5000  -120  -30]) 
ylabel (' Field  Strength  (dB)  ') 
xlabel (' Screen  Position  (m) ') 

title ([ '  Wk  =  '  num2str(Wk)  '  meters'  '  Dx  =  'num2str (Dx)  '  meters']) 


%legend ( ' Parabolic ' , ' KnifedB ' ) 

%saveas (gcf ,  ' 3D_PE_4RAY .fig'); 

%saveas (gcf,  ' 3D_PE_4RAY . bmp ' )  ; 

toe  %  Stop  stopwatch 
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